I'm looking for a fast way to generate random points in 2D according to a given 2D density function.
For instance something like this:

enter image description here

Right now I'm using a modified version of "Poisson disc" method I found here:

To make this fast, it relies on a grid on fixed interval.
In my case, I don't know the right size of the grid cells since the density changes over the plane.
(To generate this image I used no grid at all and it's pretty slow)
Is there a "right way" to generalize this algorithm for this need?
Maybe a different approach altogether?

  • $\begingroup$ How is the function specified? Is it blackbox or do you know it in advance? $\endgroup$ – Stefano Gogioso Jun 27 at 22:18
  • 1
    $\begingroup$ Do you want the samples to be independent? The fact that you're using a Poisson disc method suggests you're not. It can make a big difference. $\endgroup$ – Dan Piponi Jun 27 at 23:02
  • $\begingroup$ @StefanoGogioso function is just code I can calculate at any point $\endgroup$ – shoosh Jun 28 at 5:21
  • $\begingroup$ @DanPiponi not independent. If the density is constant, they need to be more or less constant distance from each other $\endgroup$ – shoosh Jun 28 at 5:21

Assuming you can easily compute the infimum and supremum of the density function $f$ on a given square with axis-aligned sides, you can use a quadtree to perform an approximate sampling. Let's assume that $f$ is bounded.

Work in 2D Cartesian coordinates and write $S_{l, b, \delta}$ for the axis-aligned square containing all points $(x,y)$ with $l \leq x < l+\delta$ and $b \leq y < b+\delta$, where $\delta > 0$. Assume that you start with a square domain $S_{l_0, b_0, \delta_0}$, within which you want to sample points. Fix an $\epsilon > 0$ and proceed to build a quadtree as follows, starting from the quadtree with a single node $S_{l_0, b_0, \delta_0}$ (a leaf).

For every node $S_{l, b, \delta}$ of the quadtree, compute the infimum $m_{l, b, \delta} := \inf\{f(x,y) | (x,y) \in S_{l, b, \delta} \}$ and supremum $M_{l, b, \delta} := \sup\{f(x,y) | (x,y) \in S_{l, b, \delta} \}$ of the density function $f$ on the node. These two quantities give bounds for the integral of $f$ over the node: $$ m_{l, b, \delta}\delta^2 \leq \int_{S_{l, b, \delta}} f(x,y) \,d(x,y) \leq M_{l, b, \delta} \delta^2 $$

If $(M_{l, b, \delta} - m_{l, b, \delta})\delta^2 > \epsilon$ for any leaf $S_{l, b, \delta}$ of the quadtree, proceed to subdivide the leaf into four new leaves, each with side $\delta/2$: $$ S_{l+\lambda\frac{\delta}{2}, b+\beta\frac{\delta}{2}, \frac{\delta}{2}} \hspace{5mm} \text{for} \hspace{5mm} \lambda, \beta \in \{0, 1\} $$

Because $f$ is bounded, this process is guaranteed to terminate and the maximum depth for the quadtree is $O\left(\log(\epsilon^{-1})\right)$, or more specifically it is bounded above by the smallest $d \geq 0$ such that: $$ (M_{l_0, b_0, \delta_0} - m_{l_0, b_0, \delta_0}) \delta_0^2 2^{-2d} \leq \epsilon $$

The number of nodes in the quadtree can be $O\left(\epsilon^{-1}\right)$ in specific ill-behaved cases, but it can be guaranteed to be $O\left(\log(\epsilon^{-1})\right)$ if some additional regularity requirements are imposed on $f$ (I seem to recall that $f$ being Lipschitz continuous with a fixed upper bound on the Lipschitz constant is enough).

When subdivision is terminated, define an approximate integral $I_{l, b, \delta}$ for each node as follows:

  • if $S_{l, b, \delta}$ is a leaf, set the approximate integral to $I_{l, b, \delta} := M_{l, b, \delta} \delta^2$;
  • if $S_{l, b, \delta}$ is an internal node, set $I_{l, b, \delta}$ to be the sum of the approximate integrals over its children: $$ I_{l, b, \delta} := \sum_{\lambda=0}^1 \sum_{\beta=0}^1 I_{l+\lambda\frac{\delta}{2}, b+\beta\frac{\delta}{2}, \frac{\delta}{2}} $$

Because of the termination condition for quadtree subdivision, we now have the following approximation for every node $S_{l, b, \delta}$ in the quadtree: $$ \left| I_{l, b, \delta} - \int_{S_{l, b, \delta}} f(x,y) \,d(x,y)\right| \leq \epsilon \delta^2 $$

Once the quadtree is built, $\epsilon$-approximate sampling of a point according to the density function $f$ can be performed efficiently by descending the quadtree. Starting from the root, at every internal node $S_{l, b, \delta}$, randomly choose one of the four children $S_{l+\lambda\frac{\delta}{2}, b+\beta\frac{\delta}{2}, \frac{\delta}{2}}$ with probability given by: $$ \mathbb{P}(\lambda, \beta) := \frac{I_{l+\lambda\frac{\delta}{2}, b+\beta\frac{\delta}{2}, \frac{\delta}{2}}}{I_{l, b, \delta}} $$ Once a leaf is reached, sample a uniformly distributed point within it.

This produces independent and identically distributed samples for an $\epsilon$-approximation to $f$, each sample taking $O(\log(\epsilon^{-1}))$ time and using $O(\log(\epsilon^{-1}))$ samples from a uniform probability distribution on $[0,1]$ (one sample for each step of the descent from root to leaf plus two samples for the leaf).

| cite | improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.