Placing n points in a MxM square grid I am facing an apparently well-known problem: placing $n$ points in a discrete grid so that the points are 'evenly' distributed. By evenly I mean that I would like the density of points to be nearly constant throughout the grid (avoid points concentrating in certain regions as much as possible). 
This problem seems to be quite tough (I am not a mathematician so I did not know that beforehand :). I found this paper that seemed to address a similar problem, and there they mention that a good strategy if you have a grid that is $2^k \times 2^k$ (which I can make because my grid is arbitrary) is to pick the largest empty circle and place the point at its center. So my questions are:


*

*Is that so?

*If it is, how do I guarantee that if I get several largest circles I don't accumulate the points as I am placing the points online? 

*What is a good way to implement that? I am looking at matlab but it does not seem obvious to me how to compute the largest available circle every time (I would have to sweep all the available points for placement and compute this, which will take long time). But again I am not a computer science guy either.


Any hint / help would be appreciated.
Thanks!
EDIT: to be more specific based on the comments so far, I should say that I have tried a couple of things:
1) a code I found in http://www.softimageblog.com/archives/115 posted by Anton Sherwood, which uses two relatively prime numbers to generate a pattern. This can create adequate patterns for certain combinations in R^2 but when I round the result it groups the points. I though this approach was the Halton sequences because of the prime numbers but I am not sure now.
2) I used a true Halton sequence following one of the comments, but this gave a somewhat random distribution (which is what it is supposed to do)
My problem is deterministic, as I know my grid (actually I can chose it arbitrarily with some constraints) and I also know the number of points I want to place inside, so I was more looking for a deterministic approach to place the points to achieve even distribution (as in the paper I first mentioned).
Again thanks for the help and sorry for the confusion.
 A: As mentioned already, the simplest traditional solution to the $d$-dimensional which provides quite good results in 2-dimensions is to use the Halton sequence based on the first two primes numbers (2,3), and then multiply the terms by $m$ and then round to the nearest integer.  However, as you have already found this produces fairly good - but not great - results for the discrete solution.
An alternative sequence you could use is the generalization of the Weyl / Kronecker sequence. This sequence also typically uses the first two prime numbers, however, in this case they are chosen merely because the square root of these numbers is irrational. 
It seems that you have also tried to make minor improvements by carefully choosing different prime numbers for your parameters for either the Halton or Kronecker sequences. 
However, I have recently written a detailed blog post "The unreasonable effectiveness of quasirandom sequences, on how to easily create an open-ended low discrepancy sequences in arbitrary dimensions, that is much more amenable to discretization than existing existing low discrepancy sequences, such as the Halton and Kronecker sequences.   
The section in the post called "Covering" specifically deals with your question of discretizing the low discrepancy sequences.
In the following image less red squares (which indicates a unique integer latttice point) implies a more even distribution, as each red square indicates that the cell does not contain a blue point. One can clearly see how even the $R$-sequence distributes points more evenly compared to other contemporary methods.

The solution is an additive recurrence method (modulo 1) which generalizes the 1-Dimensional problem whose solution depends on the Golden Ratio. The solution to the $d$-dimensional problem, depends on a special constant $\phi_d$, where $\phi_d$ is the unique positive root of $ x^{d+1}=x+1$.
For $d=1$,  $ \phi_1 = 1.618033989... $, which is the canonical golden ratio.
And for $d=2$, $ \phi_2 = 1.3247179572... $, which  is often called the plastic constant, and has some beautiful properties. This value was conjectured to most likely be the optimal value for a related two-dimensional problem [Hensley, 2002]. 
Jacob Rus has posted a beautiful visualization of this 2-dimensional low discrepancy sequence, which can be found here.
With this special constant in hand, the calculation of the $n$-th term is now extremely simple and fast to calculate:
$$ R: \mathbf{t}_n = \pmb{\alpha}_0 + n \pmb{\alpha}  \; (\textrm{mod} \; 1),  \quad n=1,2,3,... $$
$$ \textrm{where} \quad \pmb{\alpha} =(\frac{1}{\phi_d}, \frac{1}{\phi_d^2},\frac{1}{\phi_d^3},...\frac{1}{\phi_d^d}), $$
Of course, the reason this is called a recurrence sequence is because the above definition is equivalent to 
$$ R: \mathbf{t}_{n+1} =  \mathbf{t}_{n} +  \pmb{\alpha}  \; (\textrm{mod} \; 1) $$
In nearly all instances, the choice of $\pmb{\alpha}_0 $ does not change the key characteristics, and so for reasons of obvious simplicity, $\pmb{\alpha}_0 =\pmb{0}$ is the usual choice. However, there are some arguments, relating to symmetry, that suggest that $\pmb{\alpha}_0=\pmb{1/2}$  is a better choice.
Specifically for $d=2$, $ \phi_2 = 1.3247179572... $ and so for $\pmb{\alpha}_0= (1/2,1/2) $, 
$$\pmb{\alpha} = (0.754878, 0.56984) $$ and so the first 5 terms of the canonical 2-dimensional sequence are:


*

*(0.254878, 0.0698403)

*(0.00975533, 0.639681)

*(0.764633, 0.209521)

*(0.519511, 0.779361)

*(0.274388, 0.349201)


The Python code is
# Use Newton-Rhapson-Method
def gamma(d):
    x=1.0000
    for i in range(20):
        x = x-(pow(x,d+1)-x-1)/((d+1)*pow(x,d)-1)
    return x

d=2
n=256

# m can be any number.
# In the diagram above it is chosen to be exactly sqrt of n, 
# simply to make the visualization more intuitive s
# so that ideally each cell should contain exactly one dot.
m=16

g = gamma(d)
alpha = np.zeros(d)                 
for j in range(d):
    alpha[j] = pow(1/g,j+1) %1
z = np.zeros((n, d))    
c = (np.zeros((n,d)).astype(int)  

for i in range(n):
    z = (0.5 + alpha*(i+1)) %1
    c = (np.floor(m *z)).astype(int)
print(c)

Hope that helps!
A: Check out the Halton sequence and some other low-discrepancy sequences. 
