Drawing circles

on 14 Oct 2011 by Mukund (@muks)

How does one draw a circle? You wipe the dust off your Foley & van Dam book and look up the circle drawing algorithm, or simply use Cairo. But I had to explain to someone how to draw circles and going directly to the Bresenham method wasn't reasonable. I started with a basic implementation and we discussed ways to optimize it in #gimp. The code here was written a few days ago, and is posted for anyone who searches for it.

Circles are invisible. When we say we want to draw a circle on the screen, we want to set pixels (light up equally sized rectangular bits of the screen) at integer coordinates to represent an approximation of the circle. Let's start off by writing some code which implements a main(), handles creating an image, setting pixels in it and writing it to a Portable graymap (.pgm) file (which can be viewed with programs like Eye of GNOME, GIMP, ImageMagick's display, etc.). With that out of the way, we can work on draw_circle() alone.

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

typedef struct {
  size_t         width;
  size_t         height;
  unsigned char *data;
} Image;

static Image *
image_new (size_t width,
           size_t height)
{
  Image *image;

  image = malloc (sizeof *image);
  image->width = width;
  image->height = height;
  image->data = malloc (width * height);

  return image;
}

static void
image_free (Image *image)
{
  free (image->data);
  free (image);
}

static void
image_fill (Image         *image,
            unsigned char  value)
{
  memset (image->data, value, image->width * image->height);
}

/**
 * image_set_pixel:
 *
 * Sets a pixel passed in signed (x, y) coordinates, where (0,0) is at
 * the center of the image.
 **/
static void
image_set_pixel (Image         *image,
                 ssize_t        x,
                 ssize_t        y,
                 unsigned char  value)
{
  size_t tx, ty;
  unsigned char *p;

  tx = (image->width / 2) + x;
  ty = (image->height / 2) + y;

  p = image->data + (ty * image->width) + tx;

  *p = value;
}

static void
image_save (const Image *image,
            const char  *filename)
{
  FILE *out;

  out = fopen (filename, "wb");
  if (!out)
    return;

  fprintf (out, "P5\n");
  fprintf (out, "%zu %zu\n", image->width, image->height);
  fprintf (out, "255\n");

  fwrite (image->data, 1, image->width * image->height, out);

  fclose (out);
}

static void
draw_circle (Image *image, int radius, unsigned char value);

int
main (int argc, char *argv[])
{
  Image *image;

  image = image_new (600, 600);

  image_fill (image, 0xff);
  draw_circle (image, 200, 0);
  image_save (image, "circle.pgm");

  image_free (image);

  return 0;
}

main() basically creates a 600×600 pixels image, fills it with white, and asks draw_circle() to draw a black circle of radius 200 pixels centered in the image. Note that image_set_pixel() takes signed coordinates, where (x=0, y=0) is at the center of the image. So for an image created with image_new (400, 400), calling image_set_pixel (image, 0, 0, 0x80) will set the pixel at (200, 200) to 0x80.

So let's implement various versions of draw_circle() now. The well-known equation of a circle is \(x^2+y^2=r^2\). The variables are shown relative to each other in the figure below:

Circle equation

This type of an equation is called an implicit equation. You can substitute values for \(x\), \(y\) and \(r\) to see if the equation holds true. If it does, then the point \((x,y)\) is on the circle. So the implicit equation lets you determine if a point is on the circle or not.

Scanning a circle

In the figure above, you can see that for every line parallel to the \(y\)-axis through some \(x \in (-r,r)\), the circle intersects it in two points. If we scan \(x\) from \(-r\) to \(r\) and plot the two intersection points \((x,y_1)\) and \((x,y_2)\) in every vertical scanline, at the end of the scan we'll have drawn a circle. For this, we have to determine \(y_1\) and \(y_2\) for every \(x\). The implicit equation is no good for it, because it is only useful in checking values and not for determining them. So let's rearrange terms to find \(y\):

$$ \begin{eqnarray*} x^2+y^2 & = & r^2 \\ y^2 & = & r^2 - x^2 \\ y & = & \sqrt{r^2 - x^2} \end{eqnarray*} $$

The square-root returns a positive and negative value for \(y\), so \(y_1=y\) and \(y_2=-y\), and the points at which the circle intersects a scanline through \(x\) are \((x,y)\) and \((x,-y)\). So let's write a draw_circle() based on this:

static void
draw_circle (Image         *image,
             int            radius,
             unsigned char  value)
{
  int x, y;

  for (x = -radius; x <= radius; x++) {
    y = (int) sqrt ((double) (radius * radius) - (x * x));

    image_set_pixel (image, x, y, value);
    image_set_pixel (image, x, -y, value);
  }
}

Download circle1.c. This program generates the following image:

Circle #1

This looks like a circle, but it is discontinuous near the \(x\)-axis. At those places, for every new scanline \(x_{n+1}=x_{n}+1\), \(y\) changes by more than \(1\) (\(y_{n+1}-y_n > 1\)).. sometimes far more than \(1\). Notice that there are discontinuities between angles 0 and 45 degrees, and no discontinuities between angles 45 and 90 degrees in the drawn circle.

In a circle, all points are at the same distance from the center. In the code above, we only compute a semi-circle and mirror it by plotting at \((x,y)\) and \((x,-y)\). In the same way, we can only compute points between the angles 45 degrees and 90 degrees and mirror it 8 ways by plotting \((x,y)\), \((x,-y)\), \((-x,y)\), \((-x,-y)\), \((y,x)\), \((y,-x)\), \((-y,x)\) and \((-y,-x)\). It plots 8 arcs, each of 45 degrees. So let's modify draw_circle() to do this:

static void
draw_circle (Image         *image,
             int            radius,
             unsigned char  value)
{
  int x, y;
  int l;

  l = (int) radius * cos (M_PI / 4);

  for (x = 0; x <= l; x++) {
    y = (int) sqrt ((double) (radius * radius) - (x * x));

    image_set_pixel (image, x, y, value);
    image_set_pixel (image, x, -y, value);
    image_set_pixel (image, -x, y, value);
    image_set_pixel (image, -x, -y, value);

    image_set_pixel (image, y, x, value);
    image_set_pixel (image, y, -x, value);
    image_set_pixel (image, -y, x, value);
    image_set_pixel (image, -y, -x, value);
  }
}

Download circle2.c. This program generates the following image:

Circle #2

So we're done, right? If code to draw a circle is all we need, then we are done. On the other hand, if you want to optimize it a bit so it performs acceptably on very slow processors, keep reading.

At this point I posted a link to the code in #gimp. pippin noted that filling the disk (the area inside the circle) can be much faster than drawing the circle. The square root computation is not needed as you're merely checking every pixel's coordinates to see if it lies inside the disk or outside. This just needs testing against the implicit equation of the disk \(x^2+y^2 \leq r^2\). Note how it is very similar to the implicit equation of the circle which just checks the boundary of the disk. As the implicit equation is checked for every pixel, filling the disk can be done in a ridiculously parallel manner. Here's a modified version of draw_circle() which plots a filled disk:

static void
draw_circle (Image         *image,
             int            radius,
             unsigned char  value)
{
  int x, y;

  for (y = -radius; y <= radius; y++)
    for (x = -radius; x <= radius; x++)
      if ((x * x) + (y * y) <= (radius * radius))
        image_set_pixel (image, x, y, value);
}

Download circle3.c. This program generates the following image:

Circle #3

Ideally you'd rewrite such code to inline the image_set_pixel() as the pixels are filled in an ordered manner. A good C compiler should also eliminate redundant subexpressions and move invariants out of loops, but some embedded compilers are not good. So you'd want to add temporaries in such cases:

static void
draw_circle (Image         *image,
             int            radius,
             unsigned char  value)
{
  unsigned char *buffer;
  size_t pad;
  int r2;
  int x, y;

  buffer = image->data + (((image->height / 2) - radius) * image->width);
  pad = (image->width / 2) - radius;
  r2 = radius * radius;

  for (y = -radius; y <= radius; y++) {
    int y2;
    unsigned char *p;

    y2 = y * y;
    p = buffer + pad;

    for (x = -radius; x <= radius; x++) {
      int x2 = x * x;

      if (x2 + y2 <= r2)
        *p = value;

      p++;
    }

    buffer += image->width;
  }
}

Download circle4.c which generates the same filled disk. Then the conversation in #gimp dropped away that night. I couldn't think of a way to avoid that sqrt(), and resigned I went to bed.

The next morning I woke up to the news that Steve Jobs had died. Reading the Slashdot story, I saw a comment about rounded rectangles. The linked article mentioned how Bill Atkinson had used the arithmetic progression \(n^2=1+3+5+7+\cdots+(2n-1)\). Of course! We don't need exact square roots, only the nearest integer root. And between the angles 45 and 90 degrees where we work, \(y_{n}-y_{n+1} \leq 1\). So we can estimate whether \(y_{n+1}\) is a decrement by comparing \(y_{n+1}^2\) with \(y_n^2-(2y_n-1)\). Awesome, eh? This progression \(n^2=1+3+5+7+\cdots+(2n-1)\) was in the Knuth book, but somehow this and the circle case both didn't take their clothes off and jump into the same pool to get me all excited. Rather than give up that something can't be solved further, it's better to think that there's always a solution. It was there from before we were born just waiting to be found.

Squares progression

Also, let's see what happens to \(y=\sqrt{r^2-x^2}\) when \(x_{n+1}=x_n+1\) (the next vertical scanline).

$$ \begin{eqnarray*} y_{n+1}^2 & = & r^2 - x_{n+1}^2 \\ & = & r^2 - (x_n+1)^2 \\ & = & r^2 - (x^2_n+2x_n+1) \\ & = & r^2 - x_n^2 - 2x_n - 1 \\ & = & y_n^2 - 2x_n - 1. \end{eqnarray*} $$

So let's update draw_circle() to implement this:

static void
draw_circle (Image         *image,
             int            radius,
             unsigned char  value)
{
  int x, y;
  int l;
  int r2, y2;
  int y2_new;
  int ty;

  /* cos pi/4 = 185363 / 2^18 (approx) */
  l = (radius * 185363) >> 18;

  /* At x=0, y=radius */
  y = radius;

  r2 = y2 = y * y;
  ty = (2 * y) - 1;
  y2_new = r2 + 3;

  for (x = 0; x <= l; x++) {
    y2_new -= (2 * x) - 3;

    if ((y2 - y2_new) >= ty) {
      y2 -= ty;
      y -= 1;
      ty -= 2;
    }

    image_set_pixel (image, x, y, value);
    image_set_pixel (image, x, -y, value);
    image_set_pixel (image, -x, y, value);
    image_set_pixel (image, -x, -y, value);

    image_set_pixel (image, y, x, value);
    image_set_pixel (image, y, -x, value);
    image_set_pixel (image, -y, x, value);
    image_set_pixel (image, -y, -x, value);
  }
}

Download circle5.c which generates the same circle. Notice how only additions, subtractions and shifts are enough in the loop. image_set_pixel() can be inlined here too, but I'm not going to do it.

So how do we draw rectangles with rounded corners? We currently draw the circle as 8 arcs, right? With the appropriate border radius, we simply translate and draw 2 arcs each at the 4 corners of the rectangle. :)