10

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.

circle

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?

Ulf Gjerdingen
  • 1,414
  • 3
  • 16
  • 20
ChristopherC
  • 1,635
  • 16
  • 31
  • 5
    you should add a link to the paper you are using ... otherwise anyone want to help have to search for it and then hoping of finding the same paper and not some derivate ... discouraging most of us to actually help – Spektre Jun 03 '16 at 07:51
  • I wished but it doesn't seem to be a free paper: http://dl.acm.org/citation.cfm?id=122734&dl=ACM&coll=DL&CFID=794955035&CFTOKEN=47936980. – ChristopherC Jun 04 '16 at 05:59

2 Answers2

7

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.

daphtdazz
  • 7,754
  • 34
  • 54
  • Thanks, it's great to see such a detailed review! As said in the question, this implementation is only a naive one directly translated from the paper so I hoped that people would have referred to it while going through my code in order for it to make more sense. My own implementation actually differs from the one in the question and doesn't draw the “special case” pixel, removing the issue in this area. 1/2 – ChristopherC Jun 09 '16 at 09:45
  • As for the optimisations that you've noted, they weren't intended as such—I was concerned of having a pixel being drawn twice, possibly resulting in a wrong alpha blending, but admittedly I didn't think yet if it made sense at all or not. Anyhow, I did test your implementation and I still have the same “spike” issue in the `(sqrt(2) / 2, sqrt(2) / 2)` area. 2/2 – ChristopherC Jun 09 '16 at 09:46
  • Hmm, the last (most opaque) points that get output by my algorithm (for radius of 1000) are 705:709 @ 79.269603, 706:708 @ 78.816728, 707:707 @ 78.645375, which indicates a 45 degree gradient, which I'd expect at that location, whereas before it was a very clear horizontal line (constant j coordinate). What values are you using for radius and max opaque? – daphtdazz Jun 09 '16 at 10:03
  • And putting the `floor` back in restores that horizontal line: 706:708 @ 78.000000, 707:708 @ 78.000000, 708:708 @ 78.000000. Probably would make sense to open a chat for it if you think I can still help :-) – daphtdazz Jun 09 '16 at 10:10
  • I have tried for a max opacity of 255 and for countless of radiuses since I implemented a zoom in the viewport (this code is meant to be used within a particle simulation playground). As for the `floor`, I am “forced” to add it to your `fade_amount` variable otherwise SDL (ctypes) complains for not being passed the right type. I'm feeling bad for having already taken so much of your time but a chat would be great if you have an idea in mind! :) I'm just not too sure how to use this feature—would I need to create a new room? – ChristopherC Jun 09 '16 at 10:34
  • Quick check: the `floor` call is fine to do for the benefit of the API, but you must do it after checking whether `fade_amount` has dropped and assigning `last_fade_amount = fade_amount`; were you doing it before or after? [edit] I'll update my answer to illustrate. – daphtdazz Jun 09 '16 at 11:08
  • I am a complete idiot. Despite your detailed explanation and even the way you named your variables that make this mistake obvious, I still jumped into it. This is the proof that I didn't take the time to go through this properly. Might the 3 week trip starting tomorrow serve me as an excuse to slightly forgive my sin :') In any case, this is now working great, thanks again! Makes me wonder where I'm misunderstanding the paper—even now I'm can't tell. – ChristopherC Jun 09 '16 at 12:23
  • Note that when drawing the points as you do towards the end of the loop, this creates an almost imperceptible artifact where a couple of pixels have a higher alpha value than they should. Slightly altering these 2 line seems to fix the issue (but draws a circle 1 pixel larger), see here: http://imgur.com/IKT6WRg – ChristopherC Jun 09 '16 at 12:24
  • Now I'm left with the floating-point issue as described in an other comment: “I asked about the `float` inputs because I would like to draw a particle system and I am concerned that particles with very small velocity might look jittery if suddendly “teleporting” from one pixel to the other ever 20 frames or so”. I guess it should be another question but would you have an idea by any chance? :) It looks like I'll also need to add clipping because it's already getting very slow with a single circle and a large magnification. – ChristopherC Jun 09 '16 at 12:28
  • Addendum to the previous note: when altering the code this way, it creates a flickering pixel in some cases. Hmm. – ChristopherC Jun 09 '16 at 12:34
  • I think I answered that at the bottom of my answer... yes, I would also expect jittery teleportation with this algorithm, because it ignores the origin. I suspect that the anomalous pixels in your linked image is indeed because we are painting to the same pixel twice, once on the traversal from x=0, y=radius to x=radius/sqrt(2), y=radius/sqrt(2), and once from y=0, x=radius to that point. So you were right that you do need the guards against writing to the same locations twice! But I think at this stage I'm going to bow out :-). Don't forget to accept an answer! – daphtdazz Jun 09 '16 at 12:41
2

I got the following to work on my machine:

def draw_antialiased_circle(renderer, position, radius, alpha_max=sdl2.SDL_ALPHA_OPAQUE):

    def _draw_points(i, j):
        '''Draws 8 points, one on each octant.'''
        #Square symmetry              
        local_coord = [(i*(-1)**(k%2), j*(-1)**(k//2)) for k in range(4)]
        #Diagonal symmetry
        local_coord += [(j_, i_) for i_, j_ in local_coord]
        for i_, j_ in local_coord:
            sdl2.SDL_RenderDrawPoint(renderer, position.x + i_, position.y + j_)

    def set_alpha_color(alpha, r=255, g=255, b=255):
        '''Sets the render color.'''
        sdl2.SDL_SetRenderDrawColor(renderer, r, g, b, alpha)

    i, j, T = radius, 0, 0

    #Draw the first 4 points
    set_alpha_color(alpha_max)
    _draw_points(i, 0)

    while i > j:
        j += 1
        d = math.ceil(math.sqrt(radius**2 - j**2))

        #Decrement i only if d < T as proved by Xiaolin Wu
        i -= 1 if d < T else 0

        """
        Draw 8 points for i and i-1 keeping the total opacity constant
        and equal to :alpha_max:.
        """
        alpha = int(d/radius * alpha_max)
        for i_, alpha_ in zip((i, i-1), (alpha, alpha_max - alpha)):
            set_alpha_color(alpha_)
            _draw_points(i_, j)
        #Previous d value goes to T
        T = d

The above is a simple implementation of the diagram below.

Fig. 1

There is a mistake on this diagram: the bottom-most test should be D(r,j) < T not D(r,j) > T which may explain why you got a bit confused.

To answer your second question, I doubt you can pass some float number as a position and radius since ultimately the C function SDL_RenderDrawPoint expect int parameters.

You could pass some floats to your Python function but you will have to cast them to int before using SDL_RenderDrawPoint. There is unfortunately nothing like a half-pixel yet ;)

Jacques Gaudin
  • 15,779
  • 10
  • 54
  • 75
  • Thanks Jacques! It seems that the condition `while i > j` could be rewritten to `while i > j + 1` but in any case it's working great! Now I'm trying to understand the differences and it seems to boil down to your computation of D(r, j) differing from what I read on the paper and that I translated in my code as the variable `d`. In my version of the paper (which seems to be older than yours), I read `D(r, j) = floor((2^m - 1) * (ceil(sqrt(r^2 - j^2)) - sqrt(r^2 - j^2)) + 0.5)`, with `m = 255` here, but your version seems to only be `D(r, j) = ceil(sqrt(r^2 - j^2))`. 1/2 – ChristopherC Jun 09 '16 at 09:43
  • Your computation of the alpha value also differs and relies on the radius, which confuses me even more. If you don't mind expanding on these 2 points, I would greatly appreciate! And finally, I asked about the `float` inputs because I would like to draw a particle system and I am concerned that particles with very small velocity might look jittery if suddendly “teleporting” from one pixel to the other ever 20 frames or so. Not sure how to deal with this? 2/2 – ChristopherC Jun 09 '16 at 09:44
  • Hold on, it looks like I talked too fast. The resulting circle is without glitches but is not antialiased! – ChristopherC Jun 09 '16 at 09:54
  • Thanks for the feedback. I had a look into the alpha computation yesterday. I will elaborate on this and other details. It is actually a very good question you asked ! – Jacques Gaudin Jun 09 '16 at 10:03