A while ago a friend of mine asked me an interesting programming question: "How would you implement a DrawCircle() function?". "I'd loop and invoke Sin() and Cos() in order to obtain coordinates" I answered. "OK", he replied, "and how would you make it fast?".
This kept me thinking. So why not give it a try right here and now! Let's start with something like this (C#, .NET 2.0):
private static void DrawCircle1(int radius)
{
double angle = 0;
do
{
double rad = angle * Math.PI / 180.0;
double x = Math.Cos(rad);
double y = Math.Sin(rad);
DrawPixel(x, y);
angle++;
} while (angle < 360);
}
Sure enough, a circle appears on the screen. It looks weird - the line is dotted, not solid. Right, incrementing the angle by 1 each time might imply advancing too fast. An increment value better suited is the angle between two pixels on the circle. Something like this:
private static void DrawCircle2(int radius)
{
double angle = 0;
double step = Math.Asin(1.0 / radius) * 180.0 / Math.PI;
do
{
double rad = angle * Math.PI / 180.0;
double x = Math.Cos(rad) * radius;
double y = Math.Sin(rad) * radius;
DrawPixel(x, y);
angle += step;
} while (angle < 360);
}
This works, so that's something to start with. To get an idea about its speed, let's ignore the call to DrawPixel() for a second (and assume that we have a lightning-fast DrawPixel() implementation at that point - and by lightning-fast I mean an inlined function or macro that just does something like an indexed array lookup, some pointer arithmetic and one or two bitwise operations, not slow stuff like
SetPixel(),
Bitmap.SetPixel() or
BufferedImage.setRGB()).
Some profiling tells us that his function can calculate the pixels for 16,000 circles (of radius 100) per second.
Even with floating point units, floating point operations are expensive. Sure we can build up a large lookup table with pre-calculates Cos and Sin values, and interpolate in between lookup table values. I used to do that on other occasions, e.g. in
my 3D chat project. But here there is no way of telling in advance how precise the lookup table needs to be, and with increasing precision the table will grow in size too. It is a valid tuning measure nonetheless, but we can expect the floating point unit to provide a similar kind of optimization already wired in hardware.
Another way of speeding up the algorithm is to take advantage of the fact that a circle is symmetric. Wouldn't it be enough to calculate one quarter of the circle, and then just mirror it on the x- and y-axis. But wait - it's what about calculating only one eighth of the circle, and mirroring it on the diagonal, too? No sooner said than done:
private static void DrawCircle3(int radius)
{
double angle = 0;
double step = Math.Asin(1.0 / radius) * 180.0 / Math.PI;
do
{
double rad = angle * Math.PI / 180.0;
double x = Math.Cos(rad) * radius;
double y = Math.Sin(rad) * radius;
DrawPixel(x, y);
DrawPixel(x, -y);
DrawPixel(-x, y);
DrawPixel(-x, -y);
DrawPixel(y, x);
DrawPixel(y, -x);
DrawPixel(-y, x);
DrawPixel(-y, -x);
angle += step;
} while (angle <= 45);
}
This gives us something like the expected 8-fold speedup (120,000 circle calculations per second to be exact). Still it seems like a dead end. Runtime complexity O(N/8) still is O(N). We need a different approach.
Let's search for a way to draw a circle without the need of invoking Sin() and Cos(). What does a circle function look like? It's
x^2 + y^2 = 1
so:
y = sqrt(1 - x^2)
Let's also take advantage of the fact that we know that x is always incremented while looping. How do we know that? Because we just calculate one eighth of the circle. We are starting at coordinate (x = 0, y = radius), hence on each step will either move right or move right and down.
private static void DrawCircle4(int radius)
{
double x = 0;
double y = radius;
do
{
DrawPixel(x, y);
DrawPixel(x, -y);
DrawPixel(-x, y);
DrawPixel(-x, -y);
DrawPixel(y, x);
DrawPixel(y, -x);
DrawPixel(-y, x);
DrawPixel(-y, -x);
x++;
//cast to get rid of decimal places
y = (int)((Math.Sqrt(1.0 - Math.Pow(x / radius, 2)) * radius));
} while (x <= y);
}
100,000 circles per second. It even got slower! The calls to Sqrt() and Pow() more than eat up what has been gained by avoiding Sin() and Cos().
Still this might open a new perspective. Having a look at two neighboring pixels on the circle, what does make them neighbors? Their x- and y-coordinates differ by 1, either x does, or y, or both. Starting at a well-known pixel, e.g. (x = 0, y = radius) we could move along the circle just be deciding whether we increment x and leave y as it is, or decrement y and leave x as it is, or increment x and decrement y.
Once more the circle function:
x^2 + y^2 = 1
which equals to:
x^2 + y^2 - 1 = 0
Drawing a perfect circle on a discrete field of pixels, there is always a deviation, an error. An error e that can be expressed like this:
x^2 + y^2 - 1 = e
The pixel we draw is either inside the perfect circle or outside. When we have moved outside the perfect circle, we have to pull back to the inside, and vice versa. In each iteration, we must keep e as small as possible.
How does e change on each iteration? Let's assume we increment x and decrement y, so the next pixel is at (x + 1, y - 1):
e = x^2 + y^2 - 1
nextE = (x + 1)^2 + (y - 1)^2 - 1 = (x^2 + 2x + 1) + (y^2 - 2y + 1) - 1
The difference between nextE and e is influenced by the following factors:
dx = (x^2 + 2x + 1) - x^2 = 2x + 1
dy = (y^2 - 2y + 1) - y^2 = -2y + 1
This is how nextE looks like when x is incremented and y is decremented:
nextE = e + dx + dy
If we decide only to increment x and leave y as it is, nextE would be:
nextE = e + dx
And in case we only decrement y, nextE is:
nextE = e + dy
At this point we just have to figure out whether e stays smallest by moving along the x-axis only, or the y-axis only, or both.
Wow, no need to invoke Sqrt() or Pow(), no need for floating point functions or values, we can do it all with integers only. Note that we could even avoid the multiplications by a factor of 2 simply by writing (x + x) and (- y - y), but compilers normally replace that by a left-shift anyway. Also, we do not invoke Math.Abs(), but get rid of minus signs on our own.
private static void DrawCircle5(int radius)
{
int e = 0;
int x = 0;
int y = radius;
do
{
DrawPixel(x, y);
DrawPixel(x, -y);
DrawPixel(-x, y);
DrawPixel(-x, -y);
DrawPixel(y, x);
DrawPixel(y, -x);
DrawPixel(-y, x);
DrawPixel(-y, -x);
int dx = 2*x + 1;
int dy = -(2*y) + 1;
int e1 = e + dx;
int e2 = e + dx + dy;
int e3 = e + dy;
e1 = e1 < 0 ? -e1 : e1;
e2 = e2 < 0 ? -e2 : e2;
e3 = e3 < 0 ? -e3 : e3;
if (e1 <= e2 && e1 <= e3)
{
x++;
e += dx;
}
else if (e2 <= e1 && e2 <= e3)
{
y--;
x++;
e += dx + dy;
}
else
{
y--;
e += dy;
}
} while (x <= y);
}
1,067,000 circle calculations a second. It cannot possibly get better, can it?
Yes it can. Similar to DrawCircle4(), when we always incremented x, we can assume here that we will never enter the last else-block in DrawCircle5().
And since we are always moving to the right incrementing x, the delta of e depends entirely on whether we also move down and decrement y. The change in the error value caused by decrementing y is (-2*y + 1), so let's choose to decrement y each time e is greater or equal to y, as this will move e closer to 0 again.
All of this allows us to get rid of several statements, and some local variables too. Also, let's do the left-shifting explicitly this time instead of multiplying by 2 (as mentioned above, we may expect our compiler to take care of that, but I have experienced cases when that did not happen, e.g. when writing -2*y instead of -(2*y)).
private static void DrawCircle6(int radius)
{
int e = 0;
int x = 0;
int y = radius;
do
{
DrawPixel(x, y);
DrawPixel(x, -y);
DrawPixel(-x, y);
DrawPixel(-x, -y);
DrawPixel(y, x);
DrawPixel(y, -x);
DrawPixel(-y, x);
DrawPixel(-y, -x);
if (e < y)
{
e += (x<<1) + 1;
}
else
{
e += (x - y + 1)<<1;
y--;
}
x++;
} while (x <= y);
}
3,200,000 circles per second. It is now 200 times faster than at the beginning.
Side note: The numbers presented were measured in debug mode. In release mode DrawCircle1(), DrawCircle2() and DrawCircle3() are faster compared to debug mode by a factor around 10, DrawCircle5() and DrawCircle6() are faster by a factor around 2, while DrawCircle4() is just insignificantly faster. Anyway as you can tell from those numbers, DrawCircle5() and DrawCircle6() still clearly outperform all other variants in release mode as well.
I guess this is a good example that there is always room for improvement, and the solution that comes to mind first probably is not the quickest.
By the way, if you know an even faster algorithm (and I am sure several of them exist, maybe with optimizations only possible at assembler level), please drop me a line, and I will post a follow-up.
