Redistributing a set of uniformly distributed numbers to an arbitrarily defined shape
Here is a simple algorithm. Suppose that the shape is given by a sequence
$$ (x_0,y_0), \ldots, (x_n,y_n). $$
Define $A_i = (x_i - x_{i-1}) \cdot (y_{i-1} + y_i)$ for $1 \leq i \leq n$, and let $A = \sum_{i=1}^n A_i$, $\alpha_i = A_i/A$, and $\beta_i = \sum_{j=1}^i \alpha_i$.


  Sample a value $x$ uniformly on $[0,1]$, and determine $1 \leq i \leq n$ such that $\beta_{i-1} \leq x &lt; \beta_i$. Let $c = (x-\beta_{i-1})/(\beta_i-\beta_{i-1})$.


Note that $c$ is uniform on $[0,1]$.

Here is what we have accomplished so far:


The probability of obtaining a specific $i$ is proportional to the probability that a random point generated according to the desired density lies between $x_{i-1}$ and $x_i$.


This is because the area below the density curve between $x_{i-1}$ and $x_i$ is $(x_i - x_{i-1}) \frac{y_{i-1}+y_i}{2}$.

Consider now the function
$$
\begin{align*}
\varphi_i(t) &amp;= \int_{x_{i-1}}^t y_{i-1} + (y_i-y_{i-1}) \frac{t-x_{i-1}}{x_i-x_{i-1}}  \, dt \\ &amp;=
y_{i-1} (t-x_{i-1}) + (y_i - y_{i-1})\frac{(t-x_{i-1})^2}{2(x_i-x_{i-1})}
\end{align*}
$$
and its normalized cousin
$$
F_i(t) = \frac{\varphi_i(t)}{\varphi_i(x_i)} =
\frac{2y_{i-1}}{y_{i-1}+y_i} \frac{t-x_{i-1}}{x_i-x_{i-1}} + \frac{y_i-y_{i-1}}{y_i+y_{i+1}} \left(\frac{t-x_{i-1}}{x_i-x_{i-1}}\right)^2,
$$
which is the CDF of your density function restricted to the interval $[x_{i-1},x_i]$. 


  Output $F_i^{-1}(c)$.


Since $c \sim U(0,1)$, the resulting element has CDF $F_i$, and so has the correct distribution. It remains to explicitly invert $F_i$. Let $a = \frac{y_i-y_{i-1}}{y_i+y_{i-1}}$, $b = \frac{2y_{i-1}}{y_i+y_{i-1}}$, and $s = \frac{t-x_{i-1}}{x_i-x_{i-1}}$ (note that $0 \leq s \leq 1$). If $F(t) = c$ then $as^2 + bs = c$, and so we can calculate the output $t$ as
$$
s = \frac{-b + \sqrt{b^2+4ac}}{2a}, \qquad
t = x_{i-1} + (x_i - x_{i-1}) s.
$$
The formula for $s$ simplifies to
$$
s = \frac{\sqrt{(1-c) y_{i-1}^2 + c y_i^2} - y_{i-1}}{y_i-y_{i-1}}.
$$
We can thus rephrase the second step as


  Output
  $$ x_{i-1} + (x_i - x_{i-1}) \frac{\sqrt{(1-c) y_{i-1}^2 + c y_i^2} - y_{i-1}}{y_i-y_{i-1}}. $$




Another option for the second step is rejection sampling. Let $A$ be the area below the rectangle bounded by $(x_{i-1},0),(x_{i-1},y_{i-1}),(x_i,y_i),(x_i,0)$, and let $R$ be the upper bounding rectangle $(x_{i-1},0),(x_{i-1},y),(x_i,y),(x_i,0)$, where $y = \max(y_{i-1},y_i)$. Generate random points $(x,y)$ in $R$ repeatedly until one of them belongs to $A$, and output its ordinate $x$. Since the area of $A$ is at least half the area of $R$, on average this requires at most two samples.



In the first step you have to determine $i$ such that $\beta_{i-1} \leq x &lt; \beta_i$. While a naive implementation would take time $O(n)$, or $O(\log n)$ using binary search trees, you can use the alias method to reduce sampling time to $O(1)$.