# Optimal sampling for an arbitratry area

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.

-
Is there no interesting thing to say about the shape of the area? Do you consider that having points on the border is acceptable, or should the border be "repulsive"? Either or, without further specification, the best method is probably to draw random positions/velocities for each point inside the area, and adopt a physical model in which the points will interact together repulsively until an equilibrium is reached. – Sh3ljohn May 19 '13 at 12:57
The points can be located on borders as well - the only goal is to have the points distributed over the area as uniformly as possible. I think this is somehow related to segmentation problem where each point is a centroid of a segment. In other words, each points should have roughly same space around it so yes, the borders should be repulsive in that way. – Libor May 19 '13 at 13:04
The problem arises in a field of image processing: Given arbitrary image area, find n sample points that are as descriptive as possible about that area, while not knowing anything about the image content. – Libor May 19 '13 at 13:06
Segmentation is something else :) Regarding the points on the border, I'll ask my question differently; if there are only two points, do you want them to be at the extremity of the longest rope, or inside the surface, as far from the border AND from e/o as possible? – Sh3ljohn May 19 '13 at 13:08
Both answers are acceptable, but inside the shape is better. This is like choosing two points from interval <0,1> - I can choose either {0, 1} or {0.25, 0.75} - both answers are good since I will be working with many (hundreths or thousands) points and treating the border is negligible since only small amount of points will be located near borders. – Libor May 19 '13 at 13:55

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:

1. 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).

2. 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.

3. 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).

-