Finally I've managed to get results.
One way of generating point distributions with blue noise properties
is by means of a Poisson-Disk Distribution
Following the algorithm from the paper Fast Poisson disk sampling in
arbitrary dimensions, Robert Bridson I've got:

The steps mentioned in the paper are:
Step 0. Initialize an n-dimensional background grid for storing
samples and accelerating spatial searches. We pick the cell size to be
bounded by r/sqrt(n), so that each grid cell will contain at most one
sample, and thus the grid can be implemented as a simple n-dimensional
array of integers: the default −1 indicates no sample, a non-negative
integer gives the index of the sample located in a cell
Step 1. Select the initial sample, x0, randomly chosen uniformly from
the domain. Insert it into the background grid, and initialize the
“active list” (an array of sample indices) with this index (zero).
Step 2. While the active list is not empty, choose a random index from
it (say i). Generate up to k points chosen uniformly from the
spherical annulus between radius r and 2r around xi. For each point in
turn, check if it is within distance r of existing samples (using the
background grid to only test nearby samples). If a point is adequately
far from existing samples, emit it as the next sample and add it to
the active list. If after k attempts no such point is found, instead
remove i from the active list.
Note that for simplicity I've skipped step 0. Despite that the run-time is still reasonable. It's < .5s. Implementing this step would most definitely increase the performance.
Here's a sample code in Processing. It's a language built on top of Java so the syntax is very similar. Translating it for your purposes shouldn't be hard.
import java.util.List;
import java.util.Collections;
List<PVector> poisson_disk_sampling(int k, int r, int size)
{
List<PVector> samples = new ArrayList<PVector>();
List<PVector> active_list = new ArrayList<PVector>();
active_list.add(new PVector(random(size), random(size)));
int len;
while ((len = active_list.size()) > 0) {
// picks random index uniformly at random from the active list
int index = int(random(len));
Collections.swap(active_list, len-1, index);
PVector sample = active_list.get(len-1);
boolean found = false;
for (int i = 0; i < k; ++i) {
// generates a point uniformly at random in the sample's
// disk situated at a distance from r to 2*r
float angle = 2*PI*random(1);
float radius = random(r) + r;
PVector dv = new PVector(radius*cos(angle), radius*sin(angle));
PVector new_sample = dv.add(sample);
boolean ok = true;
for (int j = 0; j < samples.size(); ++j) {
if (dist(new_sample.x, new_sample.y,
samples.get(j).x, samples.get(j).y) <= r)
{
ok = false;
break;
}
}
if (ok) {
if (0 <= new_sample.x && new_sample.x < size &&
0 <= new_sample.y && new_sample.y < size)
{
samples.add(new_sample);
active_list.add(new_sample);
len++;
found = true;
}
}
}
if (!found)
active_list.remove(active_list.size()-1);
}
return samples;
}
List<PVector> samples;
void setup() {
int SIZE = 500;
size(500, 500);
background(255);
strokeWeight(4);
noLoop();
samples = poisson_disk_sampling(30, 10, SIZE);
}
void draw() {
for (PVector sample : samples)
point(sample.x, sample.y);
}
However, it needs to work for an infinite plane.
You control the size of the box with the parameter size
. r
controls the relative distance between the points. k
controls how many new sample should you try before rejecting the current. The paper suggests k=30
.