Fast algorithm to distribute points evenly in a 2D grid

Question

In an NxN grid, I want to select M points ($1 \leq M \leq N^2$) so that they are distributed as evenly as possible, spread out everywhere, edge to edge. Can you suggest a fast algorithm for distributing points evenly over a grid?

I won't yet try to formally define "evenly", because it's not an exact requirement. Low runtime is the most important factor here, so algorithms like "do a random distribution and refine until it satisfies a criteria" are bad, I prefer doing it in a simple nested loop.

---

Examples

Let me illustrate with some examples.

When M is a square number, the solution is simple, let's say N=5 and M=9:

trivial
#.#.#
.....
#.#.#
.....
#.#.#

It can also work in some special cases when M can be written as a sum of squares, and there's no interference between the squares, for example N=5 and M=13=9+4:

two square patterns
#.#.#
.X.X.
#.#.#
.X.X.
#.#.#

But this wouldn't exactly work for M=10, because after you select 9 points in a square pattern, there's no good place to put the 10th one.

My first shot at this was an image dithering algorithm. It looked pretty good when M was big enough, but with small M's, the error distribution makes it biased towards the edges. For example with N=3 and M=1, dithering yields a bad result when a good result is pretty obvious:

bad  good
...  ...
...  .#.
..#  ...

Let me show another example of a bad and good distribution with N=5 and M=5:

bad    good   bad    good
#....  #...#  .....  .#...
.#...  .....  .#.#.  ....#
..#..  ..#..  ..#..  ..#..
...#.  .....  .#.#.  #....
....#  #...#  .....  ...#.

I'll tell you an idea for a formal definition of "evenly". But as I said, it doesn't have to be exactly like that, it's more of an intuitive thing. So how about this: let's take every 2x2, 3x3, ... N-1xN-1 subsquare and count the point density in each of them. The deviation of those numbers should be as small as possible. Let's see for N=4 and M=4:

one   two
#.#.  #..#
....  ....
#.#.  ....
....  #..#

One:

• 2x2 sums: 1/4, 1/4, 1/4, 1/4, 1/4, 1/4, 1/4, 1/4, 1/4
• 3x3 sums: 4/9, 2/9, 2/9, 1/9
• stddev: 0.06991

Two:

• 2x2 sums: 1/4, 0, 1/4, 0, 0, 0, 1/4, 0, 1/4
• 3x3 sums: 1/9, 1/9, 1/9, 1/9
• stddev: 0.10758

(Even though number two has a higher deviation, it is spread edge to edge, so I prefer it a little, but number one works as well. Dealing with the edges is a not-so important edge-case.)

Answer 1

The optimal centroidal Voronoi tessellation for the infinite 2D grid is a hexagon tiling.

You can get one of these by placing points at even coordinates on even rows and odd coordinates for odd rows. Depends on the number of points you need.

Start with that then run Lloyd's algorithm for a few iterations to tidy up the corners and edges.

I would compute it offline on a unit square to good precision, then dilate it to your chosen rectangle with a simple linear transformation.