Rasterizing Circles and Ellipses

You might recall my post about a fast way of determining what grid squares fit within a circle of a given radius.  Beyond the roguelike applications, this algorithm also clearly useful for rasterizing a filled circle.  Here’s a modified version of the function:

def filled_circle(surface, x0, y0, r, c):
    for x in range(-r, r+1):
        dy = int(math.sqrt(r*r - x*x))
        for y in range(-dy, dy+1):
            surface.set(x0+x, y0+y, c)

surface is an object with the method set, which accepts x and y coordinates and a colour.  x0 and y0 are the coordinates of the circle’s origin, r is its radius, and c is its colour.

The function works by taking advantage of the circle equation: .  By solving for y we get the maximum value for a given x, and the function iterates between -y and y.  We can do the same thing for ellipses by substituting the ellipse equation: , where a is half the width of the ellipse and b is half the height.  Solving for y, we get the function:

def filled_ellipse(surface, x0, y0, a, b, c):
    for x in range(-a, a+1):
        dy = int(math.sqrt(b*b * (1.0 - float(x*x)/float(a*a))))
        for y in range(-dy, dy+1):
            surface.set(x0+x, y0+y, c)

That’s not a big change.

An unfilled circle uses a rather different algorithm.  If all we did was take the filled circle algorithm and for each x set pixels at dy and -dy, one would find that some pixels are missing.  This occurs in cases where there is more than one pixel for a given x.  We could account for these by swapping x and y in the function, but that results in plenty of redundant pixel setting.  Our goal is to set any given pixel only once.

A better method requires than one only finds the pixels for one-eighth of the circle and then rotate all of those pixels to find the rest of the circle.  This is simply done by going through all positive and negative combinations for the given coordinate.  For example, for (10, 1) we also have (-10, 1), (10, -1), and (-10, -1).  That’s four eighths of the circle, and the rest are determined by swapping x and y: (1, 10), (-1, 10), (1, -10), and (-1, -10).

Now the only issue is finding that eighth of a circle.  This is done in a similar way to finding the pixels in a filled circle.  We iterate through x coordinates, calculating the y coordinates.  We’ve reached an eighth when x is equal to or greater than y.  Unlike the filled circle, we will start at 1 instead of -r.  We could start at 0, but by doing that we will set the furthest top, bottom, left, and right pixels more than once.  Instead, we will set those before we start the loop.  Similarly, we will stop searching for pixels in our eighth before x is equal to y, otherwise the four more pixels will be set twice.  As such, we will set those once, when x equals y.  Thus:

def circle(surface, x0, y0, r, c):
    surface.set(x0, y0 + r, c)
    surface.set(x0, y0 - r, c)
    surface.set(x0 + r, y0, c)
    surface.set(x0 - r, y0, c)

    x = 1
    while True:
        y = int(math.sqrt(r*r - x*x))
        if x < y:
            surface.set(x0+x, y0+y, c)
            surface.set(x0-x, y0+y, c)
            surface.set(x0+x, y0-y, c)
            surface.set(x0-x, y0-y, c)
            surface.set(x0+y, y0+x, c)
            surface.set(x0-y, y0+x, c)
            surface.set(x0+y, y0-x, c)
            surface.set(x0-y, y0-x, c)
        else:
            if x == y:
                surface.set(x0+x, y0+y, c)
                surface.set(x0-x, y0+y, c)
                surface.set(x0+x, y0-y, c)
                surface.set(x0-x, y0-y, c)

            break
        x += 1

There’s a lot of lines there, but it’s mostly just pixel setting.  The actual algorithm is quite simple.

Ellipses follow a similar path, but not as simple.  Like with circles, one does not have to calculate the entire shape.  However, since ellipses are only symmetrical about two axes, we must find an entire quarter.

Recall what I mentioned about the use of the filled circle algorithm for determining the edge.  If we used that, the result is missing pixels.  Thus, when we come to the point at which pixels might be missing, we switch to iterating over the y axis.

Determining the point to switch over is not simple.  With circles it would be the halfway point, where x equals y.  With an ellipse that point is variable save for one constant: the point where the tangent is 45 degrees.  Another way of putting it is where the slope is 1 or -1.

Calculus tells us that the derivative of a function presents us with a function representing the slope of the original.  Thus we solve for y in the ellipse equation and take the derivative of the result.  I will not go into the steps involved (mostly because my calculus is rusty and I found a website to do it for me).  It suffices to say that when the slope for a given x is equal to or less than -1.0, we switch.

def ellipse(surface, x0, y0, a, b, c):
    surface.set(x0, y0+b, c)
    surface.set(x0, y0-b, c)
    surface.set(x0+a, y0, c)
    surface.set(x0-a, y0, c)

    y = 0
    for x in range(1, a):
        y = int(math.sqrt(b*b * (1.0 - float(x*x)/(float(a*a))))
        if -b * x / (a * math.sqrt(a*a - x*x)) < -1:

            break

        surface.set(x0+x, y0+y, c)
        surface.set(x0-x, y0+y, c)
        surface.set(x0+x, y0-y, c)
        surface.set(x0-x, y0-y, c)

    for dy in range(1, y+1):
        dx = int(math.sqrt(a*a*(1.0 - float(dy*dy)/float(b*b))))
        surface.set(x0+dx, y0+dy, c)
        surface.set(x0-dx, y0+dy, c)
        surface.set(x0+dx, y0-dy, c)
        surface.set(x0-dx, y0-dy, c)

That’s it and without setting any pixel more than once.