Saturday, 5 May 2012

Circles

I was recently reading @dubejf's post on generating a sin lookup table, and it reminded me of a really old technique for circle drawing, inspired by the CORDIC algorithm.

A circle rasterized on a bitmap.

Circle

In the Platonic sense, a circle is the set of all points a given distance from the origin.

As always, the Platonic ideal is unobtainable.

But for those of us fortunate enough to be stuck on the inside of a computer, we're happy to settle for an approximation instead.

Normally when I'm in a raster environment, I use the excellent Bresenham circle algorithm.  It's integer-only, and runs great on embedded devices where memory access is slow, but random access is okay.

On the other hand, if I need an anti-aliased circle and have enough silicon available, I tend to use a signed distance representation.  It's much better suited for pixel shaders and is great for feature discovery when compositing.

But this post is about a different method altogether, it approximates a circle by a polygon, and then draws that polygon instead.  It works best in vector environments, or for drawing arcs, or custom line effects, or, well, any time it's easier to draw polygons than to push pixels. Like over a remote desktop connection for example.

Anyway, lets see the algorithm, and then you can decide where it best fits in your toolbox!

Computing the maximum error.

Error term


One technique I use all the time, is to limit the error to be less than half a pixel. This way my approximation is visually indistinguishable from the ideal circle. It also makes for a great free parameter if you later need to optimize at the expense of a little quality.

Lets go ahead and calculate the error term now.

For a polygon with N vertices, let δ be the angle between each vertex:

δ = 2 π / N

Then the maximum deviation between the circle and the approximating polygon is given by:


errormax = radius . (1 - cos(δ/2) )

 Now we can go ahead and bound it from above by 0.5 pixels:

radius . (1 - cos(δ/2) ) ≤ error  ≤ 0.5 (pixels)

We can combine these equations to figure out the number of vertices our polygon needs.

N  ≥  π / cos -1(1 - 0.5 / radius)  

So now we can draw an N-sided polygon to approximate our circle, with error less than half a pixel, and it should be visually indistinguishable from the platonic ideal :

Draw a circle by using an approximating polygon.

Just to get a feel for how many vertices we're talking about, here's radius (pixels) and N (vertices) for a few select values:

radius      N   
6 8
9 10
13 12
20 15
30 18
45 22
68 26

Polygon Drawing (The Easy Way)

Well we all know one easy way to draw a polygon, lets just generate the co-ordinates for each vertex in turn, and then draw lines between them.  Something like this :

def DrawEasyPolygon(n, radius):
    prevPosition = (radius, 0)
    for i in xrange(n):
        angle = (i+1) * 2 * math.pi
        x = math.cos(angle) * radius
        y = math.sin(angle) * radius
        position = (x,y)

        DrawLine(prevPosition, position)

        prevPosition = position

But can we do better than that?  If only there was some way to get rid of all those calls to math.sin and math.cos inside the loop.

A recurrence formula for the sin function

 

Lets start with the angle sum and difference formulas:
  • sin( θ + δ ) = cos δ sin θ  + sin δ cosθ
  • sin( θ - δ ) = cos δ sin θ  - sin δ cosθ
We can add these together and rearrange to reveal:
  • sin( θ + δ )  = 2 cos δ sin θ  - sin( θ - δ )
Or, rewriting as a recurrence relation :

α = 2 cos δ
xn+1 =  α xn - xn-1

Here's what that looks like in code:

Generate a sin table, using a recurrence relation.
So now we have a little bit of trigonometry to setup, but then our entire inner loop consists of a multiply, a subtract and some register shuffling.

Incidentally this particular recurrence formula has excellent error properties. It seems you can run this recurrence forever without losing magnitude, and only a tiny predictable phase shift.  (Part of it comes from this being a so-called "symplectic" method, but it seems to perform even better than that.  I have a suspicion this particular recurrence somehow squeezes out the error during it's zero-crossing.  If I have time in the future I'd love to investigate further.)

Polygon Drawing (The Good Way)


So now we can combine our recurrence relation with our previous polygon draw:

The good way of drawing a polygon.

Some whimsical polygonal circles.

Results

Well, as you can see, we have circles!  Just for fun, I alternated drawing blue and red lines to highlight the individual line segments a little more clearly.

There's probably a little more we could do here, it looks like the maxError should maybe be lowered a little, and we may want to force N to be a multiple of four to enforce symmetry.

We'd also want to trap for very small and very large circles, and maybe tweak the polygon radius to better match the circle's radius.

Also, keep in mind our error condition, if you zoom in to this image, the error will be much bigger than a pixel and much more noticeable.

This caveat is equally applicable if you're using super-sampling.  Be sure to modify your error constraints appropriately.

But I'll leave all these extra concerns for the comments section down below.




Appendix - Circle/Polygon drawing in Python

The following two functions are hereby licensed under CC0

# From http://missingbytes.blogspot.com/2012/05/circles.html
def DrawPolygon(n,centerX,centerY,radius):
    delta=math.pi/n
    twicecosdelta=2*math.cos(delta)
    x=(radius,radius*math.cos(delta))
    y=(0,radius*math.sin(delta))

    for _ in xrange(n):
        DrawLine(x[0]+centerX,y[0]+centerY,x[1]+centerX,y[1]+centerY)
        x=(x[1],x[1]*twicecosdelta-x[0]);
        y=(y[1],y[1]*twicecosdelta-y[0]);

# From http://missingbytes.blogspot.com/2012/05/circles.html
def DrawCircle(center,radius):
    maxError=0.5 # pixels
    n=int(math.ceil(2*math.pi/math.acos(1-maxError/radius)))
    DrawPolygon(n,center[0],center[1],radius)

    No comments:

    Post a Comment