2

I need an algorithm like Bresenham's circle algorithm, but with some modifications. The algorithm must visit all pixels in the radius (so essentially a fill).

  • The algorithm must start from the center of the circle
  • It must visit all points that would normally be visited (no holes)
  • It must visit each point in the circle exactly once

One technique I came up with would first determine all pixel coordinates inside the circle by just going through the rectangle of the circle and checking with Math.Sqrt if it is inside the circle. Then it would order the pixels by distance and then visit each of them.

That would be exactly what I want, with the exception of being fast.

So my questions is: Is there a fast way to do this without fetching,ordering and then visiting each pixel?

Just for clarification I do not actually want to draw onto the image, I only want to traverse them in the described order.

Riki
  • 1,776
  • 1
  • 16
  • 31
  • Does it have to visit them in a spiral fashion, or does every next pixel have to touch the previous pixel, or will any order do as long as the first item is in the center? – Vilx- Mar 02 '15 at 11:27
  • 1
    If you're using square roots, it means you have manipulated the usual circle equation a bit. Why not use the original equation, with squarings instead of square roots? They're faster. – j_random_hacker Mar 02 '15 at 11:29
  • Once you have done that, see if you can get an expression for the difference between the value at (x, y) and at (x+1, y). If this expression is simpler than calculating the value at (x+1, y) from scratch, then you now have a cheaper way to calculate the value of the pixel to the right, given that you've already calculated the value for the current pixel. – j_random_hacker Mar 02 '15 at 11:37
  • By "The algorithm must start from the center of the circle" do you mean pixels are visited in order of distance from center? How hard is this requirement? What if pixels are same distance? – Euphoric Mar 02 '15 at 12:02
  • @Vilx & Euphoric: The next visited point should always be the one with minimal distance to the center. When there are multiple points sharing the same minimal distance it doesn't matter what point is visited first. A spiral would be ok, but its not required. Btw, yes for ordering I used just the squares, obviously one does not need sqrt. – Riki Mar 02 '15 at 12:14
  • Hmm, you know what, I reconsidered. I don't think my algorithm will work after all. :( I deleted the comments. – Vilx- Mar 02 '15 at 12:25
  • Although I don't have a solution yet, one thing that I noticed is that you could use [Radix sort](http://en.wikipedia.org/wiki/Radix_sort) here. That's very fast. – Vilx- Mar 02 '15 at 12:28
  • Another thing I noticed: pixels with the same distance come in sets of 8. Assuming that the center is at (0;0), all these pixels will have the same distance: (x;y), (-x;y), (x;-y), (-x;-y), (y;x), (-y;x), (y;-x), (-y;-x). This cuts your search space down 8 times. – Vilx- Mar 02 '15 at 12:39
  • If you must traverse many circles of different radii, you could run your algorithm (find pixels inside the circle and sort them by distance) once for a maximum radius and keep that information. Your subsequent iterations then just traverse that list until a point's distance is greater than the radius, which should be fast. – M Oehm Mar 02 '15 at 12:40
  • Since there have been several incorrect attempts already, anyone who wants to suggest an algorithm, first consider this picture: http://i62.tinypic.com/6dxopf.png. The numbers are the calculated `x*x+y*y` values of all pixels within a 10 pixel radius from the center. Values between 80 and 90 are colored magenta. Will your algorithm correctly traverse all three "85" cells, and then follow with both "89" cells? – Vilx- Mar 02 '15 at 15:13
  • @Vilx- Mine does. I even included proof that it does in my answer way before you posted this comment. – Euphoric Mar 02 '15 at 15:15
  • @Euphoric - Yes, I apologize. Your solution does work correctly. :) It does rely on sorting a bit though (PriorityQueue). – Vilx- Mar 02 '15 at 15:21

3 Answers3

3

First, we can use fact, that circle can be divided in 8 octants. So we just need to fill single octant and use simple +- coordinate change to get full circle. So if we try to fill only one octant, we need to worry only about 2 directions from center : left and left top. Also, clever use of data structures like priority queue (.NET doesn't have it, so you need to find it somewhere else) and hash map can drastically improve performance.

    /// <summary>
    /// Make sure it is structure.
    /// </summary>
    public struct Point
    {
        public int X { get; set; }
        public int Y { get; set; }

        public int DistanceSqrt()
        {
            return X * X + Y * Y;
        }
    }

    /// <summary>
    /// Points ordered by distance from center that are on "border" of the circle.
    /// </summary>
    public static PriorityQueue<Point> _pointsToAdd = new PriorityQueue<Point>();
    /// <summary>
    /// Set of pixels that were already added, so we don't visit single pixel twice. Could be replaced with 2D array of bools.
    /// </summary>
    public static HashSet<Point> _addedPoints = new HashSet<Point>();

    public static List<Point> FillCircle(int radius)
    {
        List<Point> points = new List<Point>();

        _pointsToAdd.Enqueue(new Point { X = 1, Y = 0 }, 1);
        _pointsToAdd.Enqueue(new Point { X = 1, Y = 1 }, 2);
        points.Add(new Point {X = 0, Y = 0});

        while(true)
        {
            var point = _pointsToAdd.Dequeue();
            _addedPoints.Remove(point);

            if (point.X >= radius)
                break;

            points.Add(new Point() { X = -point.X, Y = point.Y });
            points.Add(new Point() { X = point.Y, Y = point.X });
            points.Add(new Point() { X = -point.Y, Y = -point.X });
            points.Add(new Point() { X = point.X, Y = -point.Y });

            // if the pixel is on border of octant, then add it only to even half of octants
            bool isBorder = point.Y == 0 || point.X == point.Y;
            if(!isBorder)
            {
                points.Add(new Point() {X = point.X, Y = point.Y});
                points.Add(new Point() {X = -point.X, Y = -point.Y});
                points.Add(new Point() {X = -point.Y, Y = point.X});
                points.Add(new Point() {X = point.Y, Y = -point.X});
            }

            Point pointToLeft = new Point() {X = point.X + 1, Y = point.Y};
            Point pointToLeftTop = new Point() {X = point.X + 1, Y = point.Y + 1};

            if(_addedPoints.Add(pointToLeft))
            {
                // if it is first time adding this point
                _pointsToAdd.Enqueue(pointToLeft, pointToLeft.DistanceSqrt());
            }

            if(_addedPoints.Add(pointToLeftTop))
            {
                // if it is first time adding this point
                _pointsToAdd.Enqueue(pointToLeftTop, pointToLeftTop.DistanceSqrt());
            }
        }

        return points;
    }

I will leave the expansion to full list on you. Also make sure borders of the octants don't cause doubling of the points.

Ok, I couldn't handle it and did it myself. Also, to make sure it has properties you desire I did simple test :

        var points = FillCircle(50);

        bool hasDuplicates = points.Count != points.Distinct().Count();
        bool isInOrder = points.Zip(points.Skip(1), (p1, p2) => p1.DistanceSqrt() <= p2.DistanceSqrt()).All(x => x);
Euphoric
  • 12,645
  • 1
  • 30
  • 44
  • Are you sure that within the same octant all equidistant pixels will be adjacent to each other? I think that not... – Vilx- Mar 02 '15 at 14:53
  • @Vilx- What? What does "All equidistant pixels" mean? – Euphoric Mar 02 '15 at 14:56
  • Ack, sorry. OK, here's a picture. http://i62.tinypic.com/6dxopf.png I calculated the `x*x+y*y` for all pixels within a 10 pixel radius. Values between 80 and 90 have been colored magenta. Note that even within the same octant the three values of "85" are seriously apart from each other. – Vilx- Mar 02 '15 at 15:07
  • @Vilx- Yes, but that is not requirement on the algorithm. Only requirement is that points must be ordered by distance from center. There is nothing that says they need to be next to each other. – Euphoric Mar 02 '15 at 15:11
  • No, but your algorithm will fail if they're not. You only consider the adjacent pixels. – Vilx- Mar 02 '15 at 15:15
  • No, wait, I take it back. I misunderstood what it does. You're right, it works correctly. – Vilx- Mar 02 '15 at 15:17
2

I found a solution that satisfies my performance needs. It's very simple, just a offset array.

    static Point[] circleOffsets;
    static int[] radiusToMaxIndex;

    static void InitCircle(int radius)
    {
        List<Point> results = new List<Point>((radius * 2) * (radius * 2));

        for (int y = -radius; y <= radius; y++)
            for (int x = -radius; x <= radius; x++)
                results.Add(new Point(x, y));

        circleOffsets = results.OrderBy(p =>
        {
            int dx = p.X;
            int dy = p.Y;
            return dx * dx + dy * dy;
        })
        .TakeWhile(p =>
        {
            int dx = p.X;
            int dy = p.Y;
            var r = dx * dx + dy * dy;
            return r < radius * radius;
        })
        .ToArray();

        radiusToMaxIndex = new int[radius];
        for (int r = 0; r < radius; r++)
            radiusToMaxIndex[r] = FindLastIndexWithinDistance(circleOffsets, r);
    }

    static int FindLastIndexWithinDistance(Point[] offsets, int maxR)
    {
        int lastIndex = 0;

        for (int i = 0; i < offsets.Length; i++)
        {
            var p = offsets[i];
            int dx = p.X;
            int dy = p.Y;
            int r = dx * dx + dy * dy;

            if (r > maxR * maxR)
            {
                return lastIndex + 1;
            }
            lastIndex = i;
        }

        return 0;
    }

With this code you just get the index where to stop from radiusToMaxIndex, then loop through circleOffsets and visit those pixels. It will cost lot of memory like this, but you can always change the datatype of the offsets from Point to a custom one with Bytes as members.

This solution is extremely fast, fast enough for my needs. It obviously has the drawback of using some memory, but lets be honest, instantiating a System.Windows.Form uses up more memory than this...

Riki
  • 1,776
  • 1
  • 16
  • 31
1

You have already mentioned Bresenhams's circle algorithm. That is a good starting point: You could start with the centre pixel and then draw Bresenham circles of increasing size.

The problem is that the Bresenham circle algorithm will miss pixels near the diagonals in a kind of Moiré effect. In another question, I have adopted the Bresenham algorithm for drawing between an inner and outer circle. With that algorithm as base, the strategy of drawing circles in a loop works.

Because the Bresenham algorithm can place pixels only at discrete integer coordinates, the order of visiting pixels will not be strictly in order of increasing distance. But the distance will always be within one pixel of the current circle you are drawing.

An implementation is below. That's in C, but it only uses scalars, so it shouldn't be hard to adapt to C#. The setPixel is what you do to each pixel when iterating.

void xLinePos(int x1, int x2, int y)
{
    x1++;
    while (x1 <= x2) setPixel(x1++, y);
}

void yLinePos(int x, int y1, int y2)
{
    y1++;
    while (y1 <= y2) setPixel(x, y1++);
}

void xLineNeg(int x1, int x2, int y)
{
    x1--;
    while (x1 >= x2) setPixel(x1--, y);
}

void yLineNeg(int x, int y1, int y2)
{
    y1--;
    while (y1 >= y2) setPixel(x, y1--);
}

void circle2(int xc, int yc, int inner, int outer)
{
    int xo = outer;
    int xi = inner;
    int y = 0;
    int erro = 1 - xo;
    int erri = 1 - xi;

    int patch = 0;

    while (xo >= y) {         
        if (xi < y) {
            xi = y;
            patch = 1;
        }

        xLinePos(xc + xi, xc + xo, yc + y);
        yLineNeg(xc + y,  yc - xi, yc - xo);
        xLineNeg(xc - xi, xc - xo, yc - y);
        yLinePos(xc - y,  yc + xi, yc + xo);

        if (y) {
            yLinePos(xc + y,  yc + xi, yc + xo);
            xLinePos(xc + xi, xc + xo, yc - y);
            yLineNeg(xc - y,  yc - xi, yc - xo);
            xLineNeg(xc - xi, xc - xo, yc + y);
        }

        y++;

        if (erro < 0) {
            erro += 2 * y + 1;
        } else {
            xo--;
            erro += 2 * (y - xo + 1);
        }

        if (y > inner) {
            xi = y;
        } else {
            if (erri < 0) {
                erri += 2 * y + 1;
            } else {
                xi--;
                erri += 2 * (y - xi + 1);
            }
        }
    }

    if (patch) {
        y--;
        setPixel(xc + y, yc + y);
        setPixel(xc + y, yc - y);
        setPixel(xc - y, yc - y);
        setPixel(xc - y, yc + y);
    }
}

/*
 *      Scan pixels in circle in order of increasing distance
 *      from centre
 */
void scan(int xc, int yc, int r)
{
    int i;

    setPixel(xc, yc);
    for (i = 0; i < r; i++) {
        circle2(xc, yc, i, i + 1);
    }
}

This code takes care of not visiting pixels that are in two octants by skipping coincident picels on alterante octants. (Edit: There was still abug in the original code, but it's fixed now by means of the ´patch` variable.)

There's also room for improvement: The inner circle is basically the outer circle of the previous iteration, so there's no point in calculating it twice; you could keep an array of the outer points of the previous circle.

The xLinePos functions are also a bit too complicated. There are never more than two pixels drawn in that function, usually only one.

If the roughness of the search order bothers you, you can run a more exact algorithm once at the beginning of the program, where you calculate a traversing order for all circles up to a reasonable maximum radius. You can then keep that data and use it for iterating all circles with smaller radii.

Community
  • 1
  • 1
M Oehm
  • 28,726
  • 3
  • 31
  • 42