Suppose I have a probability density function defined on a region in the plane (in my case, the pdf is of the form $f(x) = \alpha\|x\|^{-\beta}$, and the region is the unit disk). For large $N$, is it possible to place $N$ points $X_1,\dots,X_N$ in the region so that the points $X_i$ are distributed according to $f(x)$, and also form a mesh of (approximately) equilateral triangles? This is clearly trivial when $f(x)$ is uniform (just put the $X_i$ in a uniform triangular lattice).

For the non-uniform case, obviously some triangles will be larger than others, but I want each individual triangle to be approximately equilateral (e.g. maximum side length and minimum side length are within 1% of each other, etc.). One possibility for the non-uniform case would be to sample $N$ points independently at random from $f(x)$ and then take their Delaunay triangulation, but I don't think there is a guarantee that the triangles will be roughly equilateral (i.e. some will be long and skinny) as $N$ becomes large.

The picture below is along these lines, if you ignore the big ugly hole in the center; each triangle is roughly equilateral, but points are not uniformly distributed.