Fitting a mesh to a density function 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.
     (source: Wayback Machine)
 A: Here is one possible interpretation of your question.

Assume a probability density function $f$ is given.
  Is there a sequence of triangulations $T_n$ with $\varepsilon_n$-equilateral triangles such that counting probability measure on nodes converges to $f$ and $\varepsilon_n\to 0$ as $n\to\infty$.

(Say a triangle is $\varepsilon_n$-equilateral if the ratio of maximum side length and minimum side length is $\le 1+\varepsilon$.)
I am almost sure that the answer is "YES" if and only if $f$ is conformal factor of a flat metric;
i.e., if and only if $f=e^{2{\cdot}\phi}$ and $\Delta \phi\equiv 0$.
A: There is an analogy in mechanics that might help: think of the nodes of the mesh as being connected by springs, which have tension proportional to something meaningful, e.g. the integral of $f(x)$ along the segment $[X_i,X_j]$. Then, if you let it stabilize, you will get a mesh with nodes distributed roughly according to $f$; if you pre-process $f$ to make is smooth enough so that it does not change much on every initial triangle, you should end up with roughly equilateral triangles, too. 
