# Fast way to generate random points in 2D according to a density function

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:

Right now I'm using a modified version of "Poisson disc" method I found here:
https://www.jasondavies.com/poisson-disc/

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?

• How is the function specified? Is it blackbox or do you know it in advance? – Stefano Gogioso Jun 27 at 22:18
• 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. – Dan Piponi Jun 27 at 23:02
• @StefanoGogioso function is just code I can calculate at any point – shoosh Jun 28 at 5:21
• @DanPiponi not independent. If the density is constant, they need to be more or less constant distance from each other – 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).