Drawing an antialiased circle as described by Xaolin Wu

Question

I'm tyring to implement the “Fast Anti-Aliased Circle Generator” routine that was described by Xiaolin Wu in his paper “An Efficient Antialiasing Technique” from Siggraph '91.

This is the code that I wrote using Python 3 and PySDL2:

def draw_antialiased_circle(renderer, position, radius):
    def _draw_point(renderer, offset, x, y):
        sdl2.SDL_RenderDrawPoint(renderer, offset.x - x, offset.y + y)
        sdl2.SDL_RenderDrawPoint(renderer, offset.x + x, offset.y + y)
        sdl2.SDL_RenderDrawPoint(renderer, offset.x - x, offset.y - y)
        sdl2.SDL_RenderDrawPoint(renderer, offset.x + x, offset.y - y)

    i = 0
    j = radius
    d = 0
    T = 0

    sdl2.SDL_SetRenderDrawColor(renderer, 255, 255, 255, sdl2.SDL_ALPHA_OPAQUE)
    _draw_point(renderer, position, i, j)

    while i < j + 1:
        i += 1
        s = math.sqrt(max(radius * radius - i * i, 0.0))
        d = math.floor(sdl2.SDL_ALPHA_OPAQUE * (math.ceil(s) - s) + 0.5)

        if d < T:
            j -= 1

        T = d

        if d > 0:
            alpha = d
            sdl2.SDL_SetRenderDrawColor(renderer, 255, 255, 255, alpha)
            _draw_point(renderer, position, i, j)
            if i != j:
                _draw_point(renderer, position, j, i)

        if (sdl2.SDL_ALPHA_OPAQUE - d) > 0:
            alpha = sdl2.SDL_ALPHA_OPAQUE - d
            sdl2.SDL_SetRenderDrawColor(renderer, 255, 255, 255, alpha)
            _draw_point(renderer, position, i, j + 1)
            if i != j + 1:
                _draw_point(renderer, position, j + 1, i)

This is a naive implementation of what I believe is being described in his paper, at the exception that I assigned the radius value to j instead of i because either I misunderstand something or there's a mistake in his paper. Indeed, he initializes i with the radius value, j with 0, and then defines the loop condition i <= j which can only be true when the radius is 0. This change led me to make some other minor modifications from what's described, and I also changed if d > T to if d < T simply because it looked broken otherwise.

This implementation works mostly well except at the start and the end of each octant, where some glitches appear.

The circle above has a radius of 1. As you can see at the start of each octant (such as in the (0, 1) area), the pixels drawn within the loop are not aligned to the first pixel that is drawn before the loop starts. Also something goes awfully wrong towards the end of each octant (such as in the (sqrt(2) / 2, sqrt(2) / 2) area). I managed to make this last issue disappear by changing the if d < T condition to if d <= T, but the same problem then shows up at the start of each octant.

Question 1: What am I doing wrong?

Question 2: Would there be any gotcha if I wanted the input position and radius to be floating points?

Answer 1

Let's deconstruct your implementation to find out what mistakes you made:

def  draw_antialiased_circle(renderer, position, radius):

First question, what does radius mean? It's impossible to draw a circle, you can only draw an annulus (a ring / donut), because unless you have some thickness you can't see it. Therefore radius is ambiguous, is it the inner radius, mid-point radius or outer radius? If you don't specify in the variable name, it will get confusing. Perhaps we can find out.

def _draw_point(renderer, offset, x, y):
    sdl2.SDL_RenderDrawPoint(renderer, offset.x - x, offset.y + y)
    sdl2.SDL_RenderDrawPoint(renderer, offset.x + x, offset.y + y)
    sdl2.SDL_RenderDrawPoint(renderer, offset.x - x, offset.y - y)
    sdl2.SDL_RenderDrawPoint(renderer, offset.x + x, offset.y - y)

OK, we're going to draw four places at once since the circle is symmetric. Why not 8 places at once?

i = 0
j = radius
d = 0
T = 0

OK, initialize i to 0 and j to radius, these must be the x and y co-ordinates. What are d and T? Not helped by undescriptive variable names. It's OK when copying scientists' algorithms to make them actually understandable with longer variable names you know!

sdl2.SDL_SetRenderDrawColor(renderer, 255, 255, 255, sdl2.SDL_ALPHA_OPAQUE)
_draw_point(renderer, position, i, j)

Huh? We're drawing a special case? Why are we doing that? Not sure. What does this mean though? It means fill in the square in (0, radius) with full opacity. So now we know what radius is, it's the outer radius of the ring, and the width of the ring apparently is one pixel. Or at least, that's what this special case tells us... let's see if that holds up in the generic code.

while i < j + 1:

We're going to loop until we reach the point on the circle where i > j, then stop. I.e. we're drawing an octant.

    i += 1

So we've already drawn all the pixels we care about at the i = 0 position, let's move onto the next one.

    s = math.sqrt(max(radius * radius - i * i, 0.0))

I.e. s is a floating point number that is the y-component of the octant at point i on the x-axis. I.e. the height above the x-axis. For some reason we have a max in there, perhaps we're worried about very small circles... but this doesn't tell us whether the point is on the outer / inner radius yet.

    d = math.floor(sdl2.SDL_ALPHA_OPAQUE * (math.ceil(s) - s) + 0.5)

Breaking it down, (math.ceil(s) - s) gives us a number between 0 and 1.0 that. This number will increase as s decreases, so as i increases, and then once it reaches 1.0 will reset to 0.0, so in a sawtooth pattern.

sdl2.SDL_ALPHA_OPAQUE * (math.ceil(s) - s) hints we are using the sawtooth input to generate a level of opaqueness, i.e. to anti-aliase the circle. The 0.5 constant addition seems unnecessary. The floor() suggests we are only interested in integer values - probably the API only takes integer values. All this shows that d ends up being an integer level between 0 and SDL_ALPHA_OPAQUE, that to begin with starts at 0, then gradually builds up as i increases, then when ceil(s) - s moves from 1 back to 0, it drops back low again too.

So what value does this take at the start? s is almost radius since i is 1, (assuming a non-trivially sized circle) so assuming we have an integral radius (which the first special case code clearly was assuming) ceil(s) - s is 0

    if d < T:
        j -= 1
    d = T

Here we spot that d has moved from being high to low, and we move down the screen so that our j position stays near to where the ring theoretically should be.

But also now we realize that the floor from the previous equation is wrong. Imagine that on one iteration d is 100.9. Then on the next iteration it has fallen, but only to 100.1. Because d and T are equal as the floor has eradicated their difference, we don't decrement j in that case, and that is crucial. I think probably explains the odd curves at the ends of your octant.

if d < 0:

Just an optimization not to draw if we're going to use an alpha value of 0

alpha = d
sdl2.SDL_SetRenderDrawColor(renderer, 255, 255, 255, alpha)
_draw_point(renderer, position, i, j)

But hang on, on the first iteration d is low, and so this is drawing to (1, radius) an almost entirely transparent element. But the special case to begin with drew an entirely opaque element to (0, radius), so clearly there's going to be a graphical glitch here. Which is what you see.

Carrying on:

if i != j:

Another small optimisation where we don't draw again if we'd be drawing to the same location

_draw_point(renderer, position, j, i)

And this explains why we only draw 4 points in _draw_point(), because you do the other symmetry here. It would simplify your code not to.

if (sdl2.SDL_ALPHA_OPAQUE - d) > 0:

Another optimisation (you shouldn't prematurely optimise your code!):

alpha = sdl2.SDL_ALPHA_OPAQUE - d
sdl2.SDL_SetRenderDrawColor(renderer, 255, 255, 255, alpha)
_draw_point(renderer, position, i, j + 1)

So on the first iteration we are writing to the position above, but at full opacity. This implies that actually it is the inner radius that we are using, not the outer radius as used by the special case error. I.e. a sort of off-by-one error.

My implementation (I didn't have the library installed, but you can see from the output checking the boundary conditions that it makes sense):

import math

def draw_antialiased_circle(outer_radius):
    def _draw_point(x, y, alpha):
        # draw the 8 symmetries.
        print('%d:%d @ %f' % (x, y, alpha))

    i = 0
    j = outer_radius
    last_fade_amount = 0
    fade_amount = 0

    MAX_OPAQUE = 100.0

    while i < j:
        height = math.sqrt(max(outer_radius * outer_radius - i * i, 0))
        fade_amount = MAX_OPAQUE * (math.ceil(height) - height)

        if fade_amount < last_fade_amount:
            # Opaqueness reset so drop down a row.
            j -= 1
        last_fade_amount = fade_amount

        # The API needs integers, so convert here now we've checked if 
        # it dropped.
        fade_amount_i = int(fade_amount)

        # We're fading out the current j row, and fading in the next one down.
        _draw_point(i, j, MAX_OPAQUE - fade_amount_i)
        _draw_point(i, j - 1, fade_amount_i)

        i += 1

Finally, would there be gotchas if you wanted to pass in a floating point origin?

Yes, there would. The algorithm assumes that the origin is at an integer / integer location, and otherwise completely ignores it. As we saw, if you pass in an integer outer_radius, the algorithm paints a 100% opaque square at the (0, outer_radius - 1) location. However, if you wanted a transform to the (0, 0.5) location, you would probably want the circle to be smoothly aliased up to a 50% opacity at both the (0, outer_radius - 1) and the (0, outer_radius) locations, which this algorithm does not give you because it ignores the origin. Therefore, if you want to use this algorithm accurately you have to round your origin before you pass it in, so there's nothing to be gained from using floats.