Fast, Antialiased Circles and Ellipses from Xiaolin Wu’s concepts

If you’ve found your way here, you may be familiar with Xiaolin Wu’s well-known line-drawing algorithm. (If not, a quick visit to wikipedia can straighten that out.)

There’s no shortage of examples and explanations for drawing straight lines using Wu’s algorithm. However, after a pretty exhaustive (or at least exhausting) search around the internet, I was unable to find any satisfactory explanation or examples for drawing circles, ellipses, or other shapes. I did find hokey programs that drew a bunch of wu-lines end-to-end to make a circle. lame.

I ended up looking up Wu’s original article in a library, and so I see it as my duty to share what I learned.

10 randomly-generated ellipses, rendered by the algorithm given below

Wu dreamt up not just an algorithm, but a principle that can be applied to drawing an anti-aliased line that moves in any varying direction. The crux of the principle is quite simple:

– if the slope is more horizontal then vertical, draw the line in a horizontal fashion, otherwise draw it in a vertical fashion.

What does it mean to draw a line in a horizontal fashion?
1. take an integer x-coordinate
2. calculate y
3. note that y should lie between two vertically adjacent pixels
4. split the full opacity between those two vertically adjacent pixels, proportional to the fractional part of y (and the inverse fraction)
5. take the next integer x-coordinate, and repeat!


for non-straight lines:
6. if the slope has changed from being more-horizontal to more-vertical, we switch to taking integer y-coords, calculating a float-valued x, and sharing the opacity between horizontally adjacent x pixels.

Once you understand the principle, you can apply it to any shape. All you have to do is calculate y from x when horizontal, and x from y when vertical. Thus, the mathematically trickiest part is often found in step 2, calculate y from an integer x (or x from an integer y). Here are a few examples:

Straight lines are very easy: y = mx + b.

Circles have the formula x^2 + y^2 = radius. Solving for y,
y = sqrt (radius^2 - x^2)

Ellipses are described by (x/a)^2 + (y/b)^2 = 1
thus,
y = b * sqrt ( 1 - (x/a)^2 )

2-d cubic splines are generally described by the following:
x(t) = a*t^3 + b*t^2 + c*t + d
y(t) = e*t^3 + f*t^2 + g*t + h

y(x) = damn near impossible

So here’s some pascal code to draw a circle in the true spirit of Wu’s principle:

First, a few types:

type
  pixel_opacity_array = array of byte;  // For use as a 2D array of opacities-- a "bitmap"!

  opacity_bitmap = record
    arr : pixel_opacity_array;
    width, height : longint;
    end;

  plot4_params_type = record
    arr : pixel_opacity_array;
    center_pix, width : longint;
    end;

const
  high_opac = high(byte);

Now for circles:

function wu_circle (const r : single {radius}) : opacity_bitmap;
  var
    xi, yi : longint; // the integer iterators. I could have used 'i' for both x and y but this was less confusing.
    xj, yj : single; // 'single' is a float type. this is the real-value y calculated from the integer x (and vice-versa)
    ffd : longint;  // the forty-five degree coordinate: stop drawing and switch from horizontal to vertical mode.
    rsq : single;  // for convenience: radius squared
    frc : single;  // the fractional part of the calculated real. its inverse will also be used.
    flr : longint;  // xj or yj rounded down (the "floor"). We'll use this twice, so give it a variable.
    p4_params : plot4_params_type;  // parameters to the actual plot function, which plots 4 points at once
  begin
    // initialization section, sets up the resultant plane and the plot4_params_type

    result.arr := nil;
    if r<=0 then exit;
    with result do begin
      width := ceil(r)*2+3;  // we need a little wiggle room in the resultant array, thus the +3
      height := width;
      setlength (arr, width**2);  // I find it easier to use dynamic arrays than 2D static arrays.
      // note, ** is the exponent operator, so width**2 = width*width.
      end;
    p.arr := result.arr;  // dynamic arrays are reference-counted pointers, no copy-on-write.
    p.width := result.width;
    with p do center_pix := (ceil(r)+1)*(width+1);
    // assuming momentarily that r is an integer, let's find center_pix (r,r)
    // first get the row: it is at r*width. Then to get the column, add r.
    // Thus (r,r) is at r*width+r = r*(width+1) -- this holds true for any r
    // r is real, so take ceil(r); also add some wiggle room, so add 1: ceil(r)+1
    // Thus: (ceil(r)+1) * (width+1)

    // Now here's the stuff you've been looking for:

    rsq := r**2;
    ffd := round (r / sqrt(2));  // forty-five-degree coord

    for xi := 0 to ffd do begin
      yj := sqrt (rsq - xi**2);  // the "step 2" formula noted above
      frc := frac(yj);
      flr := floor(yj);
      plot_4_points (xi, flr, 1-frc, p4_params);
      plot_4_points (xi, flr+1, frc, p4_params);
      end;

    for yi := 0 to ffd do begin
      xj := sqrt (rsq-yi**2);
      frc := frac (xj);
      flr := floor(xj);
      plot_4_points (flr, yi, 1-frc, p4_params);
      plot_4_points (flr+1, yi, frc, p4_params);
      end;
  end;

Here’s the plot_4_points function if you’re interested:

procedure plot_4_points (const x,y : longint; const f : single;
    var p4_params : plot4_params_type;
    const take_max : boolean = false);
  var
    opac : byte;
    yw : longint;
  procedure plot (pt : longint);
    begin
      pt := pt + p.center_pix;
      if not take_max
        then p.arr[pt] := opac
        else p.arr[pt] := max (p.arr[pt], opac);  // don't overwrite a touched pixel with a less intense value
    end;
  begin
    opac := round (high_opac*f);  // high_opac is probably 255, that is, high(byte)
    if opac=0 then exit;
    yw := y*p.width;  // I could call plot (x,y) but then I'd be calculating y*width over and over
    if (x>0) and (y>0) then begin
      plot ( x+yw);
      plot ( x-yw);
      plot (-x+yw);
      plot (-x-yw);
      end
    else if x=0 then begin
      plot ( x+yw);
      plot ( x-yw);
      end
    else if y=0 then begin
      plot ( x+yw);
      plot (-x+yw);
      end;
  end;

And finally, the ellipse routine:

function wu_ellipse (const ell_width, ell_height : single) : opacity_bitmap;
  var
    xi, yi : longint;
    xj, yj : single;
    a, b : single;  // the defining characteristics of an ellipse, width/2 and height/2
    asq, bsq : single;  // a squared and b squared
    ffd : longint;  // forty-five-degree coord
    frc : single;  // frac, as with the circle
    flr : longint;  // floor, as with the circle
    p4_params : plot4_params_type;
  begin
    // initialization stuff...

    result.arr := nil;
    if (ell_width<=0) or (ell_height<=0) then exit;
    if ell_width=ell_height then begin
      result := wucircle_to_gamma_plane(ell_width/2);
      exit;
      end;
    a := ell_width/2;
    asq := a**2;
    b := ell_height/2;
    bsq := b**2;
    with result do begin
      width := ceil(ell_width)+3;
      height := ceil(ell_height)+3;
      setlength (arr, width*height);
      p4_params.arr := arr;
      p4_params.width := width;
      p4_params.center_pix := (ceil(a)+1) + (ceil(b)+1)*width;
      end;

    // let's do this thing!

    ffd := round (asq/sqrt(bsq+asq));
    for xi := 0 to ffd do begin
      yj := b*sqrt(1-xi**2/asq);
      frc := frac(yj);
      flr := floor(yj);
      plot_4_points (xi, flr, 1-frc, p4_params);
      plot_4_points (xi, flr+1, frc, p4_params);
      end;

    ffd := round (bsq/sqrt(bsq+asq));
    for yi := 0 to ffd do begin
      xj := a*sqrt(1-yi**2/bsq);
      frc := frac (xj);
      flr := floor(xj);
      plot_4_points (flr, yi, 1-frc, p4_params);
      plot_4_points (flr+1, yi, frc, p4_params);
      end;

  end;

There are a number of things you may notice about these:

• The circle and ellipse algorithms really are very similar!
• Although the shape can have a real-valued width/height, its location must be integer. That makes it possible to plot 4 points at a time. (Actually, the circle could have plotted 8 points at a time, but then the algorithm is a little trickier to write, and more different from the ellipse one)
• (At some point I’ll implement these for real-valued locations)
• I ignored the boundary conditions! What happens at 45 degrees? Well, you can see from my screenshot that it’s just not that bad. Making the boundaries look nice can be the biggest challenge involved in using Wu’s principles, though, in some cases.
• I have not actually compiled and run the code in this particular form. Sloppy, I know. I copied and pasted it from my working project, and tweaked some things around to make it easier for y’all to read and use. Bug reports and suggestions welcome; I’ll be happy to link alternate implementations if you’ve got one.
• Thickening the line is an interesting project that will make the boundaries really ugly!

A word on thickening:
I’ve not finished my own implementation yet, but the idea is as follows. The diagram may be helpful…
– the red line indicates the thickness, which is shown perpendicular to the tangent of the line being drawn.
– thus the red line should also be approximately perpendicular to the tangents of the outer and inner boundaries
– the cyan line is the vertical strip that needs to be plotted (while moving in the more-horizontal fashion)
– the angle between them is shown in yellow

on the trigenometry of thickening a line or curve

Do a little trig, and you’ve got the length of the blue line (i.e. thickness / cosine (slope)).
I take that value and divide it by two, then alternately add and subtract that amount from the real calculated point. Then the very top and very bottom points get antialiased in Wu’s fashion.

When the line gets thick, then suddenly the boundary conditions become quite imporant! That right angle shown in green will correspond to a bite taken out of the thickened circle, where we switched from drawing horizontal lines to drawing vertical lines.

The best way I’ve thought of dealing with this is to draw lines in both directions near the boundary.