Generating points in area with at least X gap length in-between

Question

I am trying to come up with a method for generating X amount of random points in a given area (in my case a square). The one thing that makes this such an issue is that each point has to be at least Y units away from all other points.

What springs to mind at first is (in c#) to check the distance between the new point and all existing points:

while(points.Count < pointsToGenerate)
{
     Point newPoint = NewPoint();
     bool addPoint = true;
     foreach(Point p in points)
     {
          if((p - newPoint).Length() < minDistance)
          {
              addPoint = false;
              break;
          }
     }

     if(addPoint)
     {
          points.Add(newPoint);
     }
}

Now this would certainly work however if no valid points were ever found this would become a never ending loop. So throw in a magical number Z there as a limit of tries?

if(loopCount > 100)
{
     break;
}

Now this has some obvious problems. If the points are randomly generated, then loopCount can go above Z even if there are places left to place the point. It not only can, but it will happen!

What I could do is to create a list of available point per each pass and then select a random one of those. This would work flawlessly except for one thing: Performance. I'm not in need of super performance in my application, but the area is 1000^2. A lot of points to check per pass even if I limit myself to integers!

So, what I can come up with might not suffice, thus I would love some assistance on this. Is there a better way to generate X points in an area A with the minimum distance between points Y?

Thanks!

EDIT: With better, I mean generally better where a balance of performance vs perfection is achieved. A bit vague, I know. I'm not entierly sure how much overhead I can have for generating those points yet so I'm basically after anything more elegant than my own method.

~Robert

Answer 1

To understand your problem : are you looking for the optimal answer (i.e : homework), or for a very good algorithm that's better than creating random points?

In the first case, I'm afraid it's a very complicated problem if you've got no prior information on the area A. And I believe that it'll be hard to find an algorithm that's faster than exploring every single case.

However, if you've got some prior information on A then things might be simpler. For example, if it's convex, you can exploit the fact that the optimal pavement if your space is infinite is hexagonal. It means that you have to put your points (in X) at the extremities of the triangles

So :

• compute the convex hull (O(n*log(n)))
• compute the diameter (the two furthest points of your set)
• begin by adding to X one of the points (of the diameter)
• then add points that correspond to the hexagonal pavement, favoring the ones on the convex hull

This algorithm is not optimal (not unless you define a very good "favoring the ones on the convex hull"...)

---

Edit : Mr E's comment reminded me that the optimal pavement thing derives from circle packing. Kudos to him for the precision!

---

However, I do have another algorithm that looks very nice, and is perhaps even optimal! It does not require any conditions on A, and is a bit expensive, but not too much. Yeah, I know, it contradicts what I said, but who cares! Having a good algorithm is nice enough.

Let's name B the set of the available points so far. And C the points that form the extremities of B. At the beginning, B=A, and if A is a square, then C consists of 4 points (the corners). You simply have to recursively :

• compute the two furthest points of B. You only need the points in C for this
• add to X one of the points (of the diameter) at random
• remove from B the points that are now unavailable. You only need to update C for this.

I know that if you work in a 1000x1000 grid, C begins with 4 points, but after adding one point to X, it means that C grows to 1570 points (roughly (pi/2)1000). You have to notice that you never put in memory B, which is big (O(n^2) if A can be put in a nn grid) Only C, and I believe that at any point, the size of C is O(n) , which remains much better than O(n^2). And computing the diameter remains O(size(C)) = O(n)

Answer 2

Here are my thoughts, I think you have to divide the area into squares with side length equal to Y. After doing that you may be left with rectangles with one of the sides less than y, unless area = integer * y2. Now maximum number of points you can generate is number of squares+rectangles. So if X is more than that, you can end the method with failure.

To start generating the points start with the last (bottom right), smallest rectangle. Pick a random point there, and find a point on its top rectangle and left rectangle so that they are Y away from the first point, and start filling a point in a rectangle/square only if its immediate right and immediate bottom square/rectangle is filled. Hence you fill the first square at the end.

While generating the random point you only need to worry about distance from maximum two points, a point in the immediate right square/rectangle and point in the immediate bottom square/rectangle. Ofcourse if you get more than X points, you could either stop, or ignore random few points to have X points remaining

To make things more random, You can start diving the squares of side Y from any of the four sides,(here it was top left) so that the start point is different each time.

Answer 3

If limiting yourself to integers is an option, then there is an algorithm that would work. Just track the number of locations that are available in each row or column, and then recalculate those values for each row or column in the vacinity of the new point. The key here is that the random location selected is chosen only among the available points, not the entire grid.

PointsRemaining = X
Points = new CoordinateCollection
Width = 1000
Height = 1000
if (Width > Hight)
    Swap(Width, Hight)
    Swapped = true
NumberOfLocationsAvailable = new int[Width]
InitializeArrayValues(NumberOfLocationsAvailable, Height)
TotalLocationsAvailable = Width * Height
While (TotalLocationsAvailable > 0 and PointsRemaining > 0)
    NextPoint = Random(0, TotalLocationsAvailable)
    NextCoordinates = FindCoordinates(NextPoint, NumberOfLocationsAvailable)
    NumberLocationsRemoved = RemoveLocationsWithinDistance(NextCoordinates, Points, NumberOfLocationsAvailable, Y)
    TotalLocationsAvailable = TotalLocationsAvailable - NumberLocationsRemoved
    Points.Add(NextCoordinates)
    PointsRemaining = PointsRemaining - 1
if (Swapped)
    Points = RotatePoints(Points)

This puts a lot of complexity in RemoveLocationsWithinDistance, but that shouldn't be too hard. It will require an 2-dimentional square array of booleans of size Y.