I have a closed area (2D) of arbitrary shape and some number of sample points to take from the area.
Now I would like to find optimal distribution of the sample points, so that distances from sample to sample are roughly the same and the sample points cover the area as uniformly as possible:
I am doing this in a context of image processing, so the area can be defined as a set of (gray) points or by its boundary.
There may be some clever algorithm (probably iterative) to compute the positions of sample points for a given count.
Thanks for any ideas.
Here is a very simple proposal :) Note that this can totally converge to a local minimum, so if you want to be confident that you're as close as possible to the best solution, you should run the whole procedure several times, and select the one with the smallest energy... Not great, but it's not an easy problem :/
To answer your last comment, this has nothing to do with uniform sampling either: uniform sampling is randomly distributed with uniform distribution accross the shape-space. This might be suitable to initialize the current model at best, but what you're trying to achieve here is not random sampling, it's equirepartite sampling. Very much different, and extremely hard for arbitrary shapes; in fact, analytical solutions do not exist in the general case.
Dynamic Model
Let $N$ particles be defined by their 2D coordinates at time t, $p_i(t) = \big( x_i(t), y_i(t) \big)$. We define the interaction between particles $i$ and $j$ ($i\neq j$) at time $t$ as: $$ I_{i\to j}(t) = \frac{ \vec{u}_{i\to j}(t) }{ \epsilon + \| p_i(t) - p_j(t) \|_2^2 } $$
For the sake of simplicity, given a timestep $\delta t$, we suppose that changes on position and its derivatives occur only at discrete times $t_k,\ k\geq 0$. We find the following update rules for the coordinate $x$ (similarly for $y$), arbitrarily taking a unitary mass: $$ \frac{\partial x_i}{\partial t}( t_{n+1} ) = \frac{\partial x}{\partial t}( t_n ) + \delta t \sum_{j\neq i} \frac{ u_{j\to i}^x (t_n) }{ \epsilon + \| p_i(t_n) - p_j(t_n) \|_2^2 } \\ x_i(t_{n+1}) = x_i(t_n) + \delta t \frac{\partial x}{\partial t}( t_n ) + \frac{\delta t^2}{2} \sum_{j\neq i} \frac{ u_{j\to i}^x (t_n) }{ \epsilon + \| p_i(t_n) - p_j(t_n) \|_2^2 } $$
Quantitative Choices
First, the dimensions of the shape being studied should be normalized, for instance, by setting the largest dimension of the bounding box to a fixed constant $L$. $\delta t$ can be modeled as a variable of time too, recomputed at each timestep to limit the displacement of each particle to a fraction of $L$. This way, we can limit the numerical errors caused by our simplistic model.
Say we chose to limit the maximum displacement to a constant $\frac{d}{\sqrt{2}}$; from the second update equation in $x$, we find a second order polynomial in $\delta t$, which roots can be computed for each point $p_i$. Taking the smallest positive timestep found for $x$ and $y$ coordinates, we can then assign: $$ \delta t = \min( \delta t_x, \delta t_y ) $$ and finalize the update of position and velocity.
Running the Sampling
There are two additional things to consider before running the sampling:
We need a function $S$ which takes 2D coordinates and returns either; the same coordinates if they are inside the shape, OR the orthogonal projection of the input coordinates on the border of the shape otherwise. This ensures that the particles won't leave the shape. (Note: the velocity should be nulled if the coordinates are projected).
The particles should be initialized with random positions inside the shape (rejection sampling in the bounding box with the previous function $S$ should do), and random velocities (direction and magnitude). That's 5 random values per particle.
Repeat the update until the maximum displacement between two timesteps goes below a certain fixed threshold, that's you convergence criterion.
Final note: I'm absolutely not sure whether this yields accurate results for fancy shapes. If you get unsatisfactory/weird results, a straightforward extension would be to cancel the interactions between particles that are not "facing" each-other (ie, $[p_i,p_j]$ crosses the boundary).