In a 2D pixel array, I need an efficient algorithm that will select p% of pixels that are the most spread out.
This can be done adaptively by selecting points, then repeatedly adjusting the positions of points that are too close together. But this isn't efficient since it would require many iterations and distance calculations.
It doesn't have to be perfect, it just needs to avoid point clusters as much as can be done efficiently.
You want a Poisson Disk distribution, but it's tricky. Doing a search turns up lots of academic papers about how to do it efficiently: http://people.csail.mit.edu/thouis/JonesPoissonPreprint.pdf
Thanks to everyone for the answers!
The best solution appears to be using "prebuilt building blocks": n x n arrays with already-selected cells, and covering the pixel array with these.
For example, a 4 x 4 array with 12.5% coverage would be:
0 0 1 0
0 0 0 0
1 0 0 0
0 0 0 0
With 6.3% coverage:
0 0 0 0
0 1 0 0
0 0 0 0
0 0 0 0
To get a % coverage between these, just alternate between these blocks according to a running tally of the overall actual % coverage so far. To cover a width that's not a multiple of 4, use some 3 x 3 blocks. To cover a larger area more efficiently, just use larger blocks.
This covers the whole array efficiently with no distance calculations or floating-point arithmetic.
The "most spread out" selection of pixels is the set whose Delaunay triangulation consists of equilateral triangles. The set of points which leads to this triangulation is found by splitting the pixel array into a set of boxes, where each box is sqrt(3) longer than it is wide. Each box contributes 5 pixels to the final pixel set (one at each corner, plus a center node at the center of the box). The trick is to find how many rows and columns of boxes will give you this 1:sqrt(3) ratio. Without going through the derivation, here's how you get that:
std::vector<PixelIndices> PickPixels(int width, int height, float percent)
{
int total_pixels = width*height;
int desired_pixel_count = (int)total_pixels*percent;
// split the region up into "boxes" with 4 corner nodes and a center node.
// each box is sqrt(3) times taller than it is wide.
// calculate how many columns of boxes
float a = 1.155*height/(float)width;
float b = .577*height/(float)width + 1;
float c = 1 - desired_pixel_count;
int num_columns = (int)((-b + sqrt(b*b -4*a*c))/(2*a));
// Now calculate how many rows
int num_rows = .577*height*num_columns/(float)width;
// total number of pixels
int actual_pixel_count = 2*num_rows*num_columns + num_rows + num_columns + 1;
std::cout << " Total pixels: " << total_pixels << std::endl;
std::cout << " Percent: " << percent << std::endl;
std::cout << "Desired pixels: " << desired_pixel_count << std::endl;
std::cout << " Actual pixels: " << actual_pixel_count << std::endl;
std::cout << " Number Rows: " << num_rows << std::endl;
std::cout << "Number Columns: " << num_columns << std::endl;
// Pre-allocate space for the pixels
std::vector<PixelIndices> results;
results.reserve(actual_pixel_count);
// Now get the pixels, my integer math is probably wrong here, didn't test
// (didn't even finish, ran out of time)
for (int row = 0; row <= num_rows; row++)
{
int row_index = row*height/num_rows;
// Top of box
for (int col = 0; col <= num_columns; col++)
{
int col_index = col*width/num_columns;
results.push_back(PixelIndices(row_index, col_index));
}
// Middle of box
if (row != num_columns)
{
for (int col = 0; col < num_columns; col++)
{
// I'll leave it to you to get this, I gotta go!
}
}
}
return results;
}
Instead of using integer division to find the indices, you could speed this up by finding the distance between each point in a row/column and just adding by the offset.
You might try Wang tiles:
http://en.wikipedia.org/wiki/Wang_tile
(See the pages linked there to Cohen's paper, and Kopf's paper. I'm a new user so can't post all the links).
These tie together both the prebuilt tiles idea, as well as the evenly distributed requirement usually solved with Poisson-disk patterns. Wang tiles can avoid periodicity effects that are almost surely an issue with more direct use of prebuilt tiles.
Quite old, but worth a dig, because answers missed an important method and focused on optimal solutions, which you are not interested in. But the method I suggest may or may not suit your needs however.
You can use quasi random sequences, which are designed for such problems. Most widespread are Sobol sequences, for which you can find canned packages for virtually any language. They are blazingly fast: only bitwise arithmetic.
They most likely will produce some clusters, but this can be avoided by selecting the "seed" to use for the x and y dimensions beforehand, and checking with naked eye.
It depends on what you want to do with the points: if "visual spread out" is important, this may not be what you want. If you want points which "fill the plane" almost evenly, they perfectly do the job. They are especially useful to average something on an image quickly, since it requires less points that with "normal" random generation. Experiment with different dimensions and see.
See also this link for exemples experiments, and pictures.
How about this:
This may be accurate enough, but if not, you could always replace step 3 with:
"Remove the point that has the smallest sum, and if you need to remove more points to reach your desired percentage, then go back to step 1."
Wait. Now I'm wondering. Are you trying to find the points that are most spread out from a given set of points...or trying, from a given array, to find the points that would be the most spread out? That's totally different...and still really hard.
How about calculating a "density" value for each pixel to start with, based on its proximity to all the other pixels. Then, repeatedly remove the most "dense" pixel until you're below p% remaining in the list.
You'd need to do the distance calculation to determine density between any given two points at most twice. First time would be when you build the original list - each pixel would need to be compared with each other pixel. Second would be when you remove a pixel from the list - you'd have to calculate the removed one against each one remaining in the list. This is to account for density values changing as each pixel is removed - for example, 2 pixels directly next to each other would have a very very high value, but once one is removed, the remaining one might have a very low value.
Some quick pseudocode (note that in this example, higher density areas have a low number)
For I = 0 to MaxPixel
For J = I to MaxPixel
PixelDensity[I], PixelDensity[J] += DistanceBetween(Pixels[I], Pixels[J])
While PixelDensity.Count > TargetCount
RemovePixel = IndexOfSmallest(PixelDensity)
ForEach I in PixelDensity
PixelDensity[I] -= DistanceBetween(Pixels[I], Pixels[RemovePixel])
PixelDensity.Remove(RemovePixel)
If memory is less of a concern than calculation time, you can also store the distance between any two points in a simple 2d array. Also, instead of simple distance, it might be helpful to make the distance calculation exponential - that would avoid something like having two points almost on top of each other, but far away from everything else, and having both make it through.
An iterative mine-sweeper flood-fill approach would be straightforward to visualise.
Ooh! How about this!
(In a very hand-wavey way, since I don't know whether your matrix is square or whatever...I'll assume it is.)
Say you have a 1000x1000 array that you want to place 47 point into (I'm picking 47 so that it is an unusual number that won't fit "nicely").
You take the ceil(sqrt(47))...that will give you a value (7). So we make a 7x7 square, fill it with 47 pixels (some are blank), and imagine placing that in the corner of the array.
Now, translate each of those pixels to a new location, based on where they are in the small (7x7) array to the large array (1000x1000). A simple equation should do this for you...for the X coordinate, for example:
xBigArrayIndex = xSmallArrayIndex * 1000 / 7;
Then your pixels will be super spread out! And it's nice and fast.
The only downside is that this only works perfectly when your square is ideally spaced to begin with...if you fill it naively (starting at the top left, going across, etc.) you will end up with with a slightly sub-ideal spread...since the translated pixels won't quite reach the bottom right of the large array. But maybe this is good enough? And if not, perhaps it's a smaller subset of the problem that is easier to deal with?
You can Use convex hull algorithm, and exclude points which this algorithm will compute and repeat it as long as it get to your p% criteria, or
execute convex hull algorithm steps, check for points included in hull and inside of it to meet the criteria 100% - p%
some demos of convex hull are here http://www.cs.unc.edu/~snoeyink/demos/
and here you got some more info http://softsurfer.com/Archive/algorithm_0109/algorithm_0109.htm
