Midpoint circle algorithm for filled circles

Question

The Midpoint circle algorithm can be used rasterize the border of a circle. However, I want the circle to be filled, without drawing pixels multiple times (this is very important).

This answer provides a modification of the algorithm that yields a filled circle, but some pixels are visited several times: fast algorithm for drawing filled circles?

Q: How can I rasterize a circle without drawing pixels multiple times? Note that RAM is very limited!

Update:

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace CircleTest
{
    class Program
    {
        static void Main(string[] args)
        {
            byte[,] buffer = new byte[50, 50];
            circle(buffer, 25, 25, 20);

            for (int y = 0; y < 50; ++y)
            {
                for (int x = 0; x < 50; ++x)
                    Console.Write(buffer[y, x].ToString());

                Console.WriteLine();
            }
        }

        // 'cx' and 'cy' denote the offset of the circle center from the origin.
        static void circle(byte[,] buffer, int cx, int cy, int radius)
        {
            int error = -radius;
            int x = radius;
            int y = 0;

            // The following while loop may altered to 'while (x > y)' for a
            // performance benefit, as long as a call to 'plot4points' follows
            // the body of the loop. This allows for the elimination of the
            // '(x != y)' test in 'plot8points', providing a further benefit.
            //
            // For the sake of clarity, this is not shown here.
            while (x >= y)
            {
                plot8points(buffer, cx, cy, x, y);

                error += y;
                ++y;
                error += y;

                // The following test may be implemented in assembly language in
                // most machines by testing the carry flag after adding 'y' to
                // the value of 'error' in the previous step, since 'error'
                // nominally has a negative value.
                if (error >= 0)
                {
                    error -= x;
                    --x;
                    error -= x;
                }
            }
        }

        static void plot8points(byte[,] buffer, int cx, int cy, int x, int y)
        {
            plot4points(buffer, cx, cy, x, y);
            if (x != y) plot4points(buffer, cx, cy, y, x);
        }

        // The '(x != 0 && y != 0)' test in the last line of this function
        // may be omitted for a performance benefit if the radius of the
        // circle is known to be non-zero.
        static void plot4points(byte[,] buffer, int cx, int cy, int x, int y)
        {
#if false // Outlined circle are indeed plotted correctly!
            setPixel(buffer, cx + x, cy + y);
            if (x != 0) setPixel(buffer, cx - x, cy + y);
            if (y != 0) setPixel(buffer, cx + x, cy - y);
            if (x != 0 && y != 0) setPixel(buffer, cx - x, cy - y);
#else // But the filled version plots some pixels multiple times...
            horizontalLine(buffer, cx - x, cy + y, cx + x);
            //if (x != 0) setPixel(buffer, cx - x, cy + y);
            //if (y != 0) setPixel(buffer, cx + x, cy - y);
            //if (x != 0 && y != 0) setPixel(buffer, cx - x, cy - y);
#endif
        }

        static void setPixel(byte[,] buffer, int x, int y)
        {
            buffer[y, x]++;
        }

        static void horizontalLine(byte[,] buffer, int x0, int y0, int x1)
        {
            for (int x = x0; x <= x1; ++x)
                setPixel(buffer, x, y0);
        }
    }
}

Here's the relevant result:

00000111111111111111111111111111111111111111110000
00000111111111111111111111111111111111111111110000
00000111111111111111111111111111111111111111110000
00000111111111111111111111111111111111111111110000
00000111111111111111111111111111111111111111110000
00000011111111111111111111111111111111111111100000
00000011111111111111111111111111111111111111100000
00000011111111111111111111111111111111111111100000
00000001111111111111111111111111111111111111000000
00000001111111111111111111111111111111111111000000
00000000111111111111111111111111111111111110000000
00000000111111111111111111111111111111111110000000
00000000011111111111111111111111111111111100000000
00000000001111111111111111111111111111111000000000
00000000000111111111111111111111111111110000000000
00000000000011111111111111111111111111100000000000
00000000000001111111111111111111111111000000000000
00000000000000122222222222222222222210000000000000
00000000000000001222222222222222221000000000000000
00000000000000000012333333333332100000000000000000
00000000000000000000012345432100000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000

The bottom pixels are plotted too many times. What am I missing here?

Update #2: This solution works:

static void circle(byte[,] buffer, int cx, int cy, int radius)
{
    int error = -radius;
    int x = radius;
    int y = 0;

    while (x >= y)
    {
        int lastY = y;

        error += y;
        ++y;
        error += y;

        plot4points(buffer, cx, cy, x, lastY);

        if (error >= 0)
        {
            if (x != lastY)
                plot4points(buffer, cx, cy, lastY, x);

            error -= x;
            --x;
            error -= x;
        }
    }
}

static void plot4points(byte[,] buffer, int cx, int cy, int x, int y)
{
    horizontalLine(buffer, cx - x, cy + y, cx + x);
    if (y != 0)
        horizontalLine(buffer, cx - x, cy - y, cx + x);
}    

Solution 1

The answer to the other question is perfectly fine. However since it is creating confusion, I'm going to explain it a bit.

The algorithm you see in Wikipedia basically finds x and y of 1/8 of a circle (angles 0 to pi/4) and then draws 8 points which are its mirrors. For example:

    (o-y,o+x) x         x (o+y,o+x)

(o-x,o+y) x                  x (o+x,o+y) <-- compute x,y

                   o

(o-x,o-y) x                  x (o+x,o-y)

    (o-y,o-x) x         x (o+y,o-x)

What the other solution suggests, which makes perfect sense if you look closely to this picture, is to instead of drawing 8 points, draw 4 horizontal lines:

    (o-y,o+x) x---------x (o+y,o+x)

(o-x,o+y) x-----------------x (o+x,o+y) <-- compute x,y

                   o

(o-x,o-y) x-----------------x (o+x,o-y)

    (o-y,o-x) x---------x (o+y,o-x)

Now if you compute (x,y) for angles in [0, pi/4] and draw these 4 lines for every computed point, you will have drawn many horizontal lines filling a circle without any line overlapping the other.

Update

The reason you get overlapping lines in the bottom of the circle is that the (x,y) coordinates are rounded, so in those locations the (x,y) move horizontally themselves.

If you take a look at this wikipedia picture:

You will notice that on the top of the circle, some pixels are horizontally aligned. Drawing horizontal lines originating from those points overlap.

If you don't want this, the solution is quite easy. You have to keep the previous x you have drawn with (since the top and bottom are mirrors of the original (x,y), you should keep the previous x which represents the y of those lines) and only draw the horizontal lines if that value changes. If it doesn't, it means that you are on the same line.

Given the fact that you will first encounter the innermost points, you should draw lines for the previous point only if the new point has different x (of course, the last line is drawn always). Alternatively, you can start drawing from angle PI/4 down to 0 instead of 0 to PI/4 and that you will first encounter the outer points, therefore you draw lines every time you see a new x.

Solution 2

I needed to do this, here's the code I came up with for it. The visual image here shows the pixels drawn where the number is the order in which the pixels are traversed, and the green numbers represent pixels that are drawn using the reflection of the completion of a column using symmetry as shown in the code.

void drawFilledMidpointCircleSinglePixelVisit( int centerX, int centerY, int radius )
{
    int x = radius;
    int y = 0;
    int radiusError = 1 - x;

    while (x >= y)  // iterate to the circle diagonal
    {

        // use symmetry to draw the two horizontal lines at this Y with a special case to draw
        // only one line at the centerY where y == 0
        int startX = -x + centerX;
        int endX = x + centerX;
        drawHorizontalLine( startX, endX, y + centerY );
        if (y != 0)
        {
            drawHorizontalLine( startX, endX, -y + centerY );
        }

        // move Y one line
        y++;

        // calculate or maintain new x
        if (radiusError<0)
        {
            radiusError += 2 * y + 1;
        }
        else 
        {
            // we're about to move x over one, this means we completed a column of X values, use
            // symmetry to draw those complete columns as horizontal lines at the top and bottom of the circle
            // beyond the diagonal of the main loop
            if (x >= y)
            {
                startX = -y + 1 + centerX;
                endX = y - 1 + centerX;
                drawHorizontalLine( startX, endX,  x + centerY );
                drawHorizontalLine( startX, endX, -x + centerY );
            }
            x--;
            radiusError += 2 * (y - x + 1);
        }

    }

}

Solution 3

I came up with an algorithm that draws the circle already filled.
It iterates over the pixels that the circle will be drawn upon and nothing else.
From here on it all about the speed of the draw-pixel function.

Here's a *.gif that demonstrates what the algorithm does !

As for the algorithm here's the code :

    //The center of the circle and its radius.
    int x = 100;
    int y = 100;
    int r = 50;
    //This here is sin(45) but i just hard-coded it.
    float sinus = 0.70710678118;
    //This is the distance on the axis from sin(90) to sin(45). 
    int range = r/(2*sinus);
    for(int i = r ; i >= range ; --i)
    {
        int j = sqrt(r*r - i*i);
        for(int k = -j ; k <= j ; k++)
        {
            //We draw all the 4 sides at the same time.
            PutPixel(x-k,y+i);
            PutPixel(x-k,y-i);
            PutPixel(x+i,y+k);
            PutPixel(x-i,y-k);
        }
    }
    //To fill the circle we draw the circumscribed square.
    range = r*sinus;
    for(int i = x - range + 1 ; i < x + range ; i++)
    {
        for(int j = y - range + 1 ; j < y + range ; j++)
        {
            PutPixel(i,j);
        }
    }

Hope this helps ... some new users ... sorry for necro-posting.
~Shmiggy

Solution 4

I wanted to comment on your Update #2: This solution works: (but I guess I need more reputation first...) that there's a small bug in the solution, coincidentally when drawing small circles. If you set the radius to 1 you get

00000
00000
01110
00100
00000

To fix this all you need to do is change the conditional check in plot4points from

if (x != 0 && y != 0)

to

if (y != 0)

I've tested this on small and big circles to make sure each pixel is still assigned only once. Seems to work great. Me thinks the x != 0 was not needed. Save a tiny bit of performance too.

View more solutions

Comments

[Shahbaz]

The answer does not visit any pixels several times. Why do you say that?

[l33t]

My implementation keeps doing multi-plotting at the top/bottom of the circle. Maybe I just didn't understand the answer?

[l33t]

Thanks. I'll update my question with sample code that demonstrates the issue. There must be a certain condition that I'm missing.

[l33t]

Hm. So I need two "previous line" variables? One for 'y' and one for 'x' (mirrored)?

[l33t]

Well, it's not that simple. To draw a line you need the outer pixel. And the algorithm does not guarantee that the most outer pixel is visited first. So if I track visited lines, some pixels will be missed.

[Shahbaz]

@NOPslider, you need only 1, which is for x. Let's consider the same image from Wikipedia. You start at angle 0, so you have (x,0) and you start going up. As you go up, for a few pixels x remains the same and y changes. 2 of the horizontal lines will be drawn anyway and the other two (mirrored over line y = x) must be checked whether they introduce a horizontal line. If they do, you draw them.

[Shahbaz]

Also, you are right, you visit the node from inside out. I'll update the answer

[l33t]

Thanks. So I need to edit both the plot4points() function and the outer loop?

[l33t]

Ok. I obviously need to modify the main loop. Any ideas on how to reverse the "rotation" of this optimized version of Bresenham? I'm getting all dizzy by this algorithm :)

[Shahbaz]

In the first method, you need to keep last_x and last_y and draw lines (last_x, last_y) -- (-last_x, last_y) and (last_x, -last_y) -- (-last_x, -last_y) always and lines (last_y, last_x) -- (-last_y, last_x) and (last_y, -last_x) -- (-last_y, -last_x) only if last_x != x. So you need to also divide the plot4points (or better, make it draw two lines (instead of 4) and give the mirrored x and y in the second call (which is conditioned)).

[Shahbaz]

In the second method, instead of starting from (x=radius,y=0) while x >= y, you need to start from (radius*sqrt(2)/2, radius*sqrt(2)/2) while x >= 0. To avoid floating point computation, you may want to represent sqrt(2) by an equivalent a/b number, preferably with a b as a power of 2 so it could be replaced with a shift. You then need to reverse all operations. That is all -s to +s and vice versa.

[l33t]

You're being really helpful. Thanks! How many points do you want to modify my sample code? When I apply your changes, I either get weird x/y values or the drawn figure is... something else than a circle :P. I'll update my question with the proposed solution.

[Shahbaz]

I saw it. Well, you have also negated the operations on error, but not the if on error. You should either keep error as it was, or change the if to if (error <= 0). I'd say keep error is it has originally been. Honestly though, I'm not sure what the initial value of error should be in this case.

[l33t]

How about running the original while loop without plotting? Then you get the ending x=17, y=19 and error=0 and could run your reversed loop...

[l33t]

let us continue this discussion in chat

[l33t]

See my update. I solved it using a quite nice trick. I will accept your answer as you pointed me in the right direction. Thanks!

[l33t]

Thanks, but please be aware that the sqrt() call is very expensive. The same goes for multiplication and division.

[Aaron]

On x86 procs multiplication, division and sqrt are all single instructions.

[l33t]

This makes sense! And it could actually imply a noticeable performance gain. Branching is usually crazy expensive on GPUs.

[l33t]

A single instruction does not necessarily mean high performance. There are clock cycles to consider, precision issues, platform differences... It makes no sense to use mathematically complex algorithms (e.g. square root approximation) when there is no need to.