Random Points in a Triangle (Generating Random Sequences II)

In the first post of this series, I discussed a method to generate a random (guaranteed) permutation. On another occasion, I presented a method to recycle expensive random bits. Today, I will discuss something I had to do recently: generate random points inside a triangle.

Let us start by a rectangle. If I want to generate cloud of random dots in a rectangle, you would be easily convinced that generating pairs of (pseudo) random numbers, one for each coordinate, is quite sufficient. Let \mathcal{U}(a,b) be a (pseudo) random function that returns a (pseudo) random number such that a\leqslant\mathcal{U}(a,b)\leqslant{}b. So for a rectangle of height h and width w, a cloud of points generated by (\mathcal{U}(0,w),\mathcal{U}(0,h)) would look quite (uniform) random.

But what if we have, say, a triangle? It’s not quite easy, but it shouldn’t be too hard. At first, we could think that it suffices to compute the triangle’s bounding box and use the rejection method: we draw random points inside the rectangle and if a random point falls outside the triangle (while being inside the rectangle) it is rejected and we try again. This method is good if the triangle to rectangle surface ratio is good, but disastrous if the triangle is very small compared to its bounding rectangle. Say, something like:

So this first solution is not very good. A much better solution asks for some linear algebra, but not much, be assured. Let us say we have a triangle like the following:

The two arrows, or vectors, can be determined easily: take any corner as the common base, and draw them to opposite corners (in fact any two vectors in any order will do, but for clarity here, let us consider that they form a corner as in the above picture). These two vectors form, if they are different, a vector space basis. The space here is in fact only a plane, the plane in which the triangle lies. Let u and v be those two vectors.

Now, let’s generate random points in the triangle. The two vectors u and v generate the parallelogram:

but only the darker gray region is our triangle; the pale region is to be ‘rejected’. Now: let u=(u_1,u_2) and v=(v_1,v_2) be the vectors and their components. Let also (x,y)=(\mathcal{U}(0,1),\mathcal{U}(0,1) be a random point in a unit square. This point can be mapped in the parallelogram using the following transform:

\left[\begin{matrix}u_1 &v_1\\u_2&v2\end{matrix}\right]\left[\begin{matrix}x\\y\end{matrix}\right]=\left[\begin{matrix}p_1\\p_2\end{matrix}\right]

The parallelogram is therefore a sheared and scaled version of the unit square and (p1,p2) are the x,y coordinates of the random point within the parallelogram. But we want only inside the original triangle!

We can now use rejection in the original unit square to make sure we only yield points in the triangle: we can reject all pairs with x+y>1, since they will be mapped in the second (pale gray) triangle. We then draw random pairs until we get x+y\leqslant{}1, and then use the above transform to map it into the triangle.

It doesn’t take long to get a Mathematica implementation running. (I suppose you could do the same with C and Gnuplot, only it’s not as immediate.) My implementation for two random vector u and v yields the following cloud of points:

And we see that the points are uniformely distributed, showing (approximately) equal density everywhere in the triangle.

15 Responses to Random Points in a Triangle (Generating Random Sequences II)

  1. Question says:

    Couldn’t you use the fact that the two triangles of the parallelogram are mirror images to map points that fall in the x+y>1 triangle back onto the original triangle?

    • Dude says:

      Not only that, but uniformity should be preserved across affine transforms in the plane (i.e. if you can create a uniform distribution on a unit square you can transform it to any parallelogram and maintain the uniformity), so you should just generate x, y = rand(), rand() in [0..1], project onto the parallelogram and then mirror if the point is on the other side of the diagonal (no “inside” checking needed).

  2. dff says:

    What is implicit but important is that the density of the points remains uniform after the transformation by virtue of the fact that it is linear.

  3. ugasoft says:

    I’ve used random sampling over baricentric coordinates.

  4. munificent says:

    You can do that without having to ditch half your points by remapping ones outside the triangle back into it by mirroring them along the diagonal:

    vec pointInTriangle(vec a, vec b, vec c)
        // make basis vectors from a->b and a->c
        vec u = b - a
        vec v = c - a
        // pick a random point in the unit square
        vec unit = vec(rand(0, 1), rand(0, 1))
        // if the point is outside the triangle, remap it inside
        // by mirroring it along the diagonal
        if unit.x + unit.y > 1 then
            unit = vec(1 - unit.y, 1 - unit.x)
        // now transform it to fit the basis
        return vec(unit.x * u.x + unit.y * v.x, unit.y * u.y + unit.y * v.y)
  5. […] This post was mentioned on Twitter by Proggit Articles, Nazar Vinnichuk. Nazar Vinnichuk said: Random Points in a Triangle (Generating Random Sequences II) « Harder, Better, Faster, Stronger http://bit.ly/dzDka0 […]

  6. Steven Pigeon says:

    Many of you seemed to have had the same idea, and that’s a good optimization, especially as the cost of generating a random pair grows.

  7. andrei says:

    why discard half of the pairs, when instead of rand(1)+rand(1) you could have rand(0.5)+rand(0.5) and that way you’re sure you always fall below zero?

  8. […] method and I explained (although indirectly) how a linear congruential generator works. In a second post, I explained how to generate 2D points uniformly and randomly distributed a triangle, discussing […]

  9. […] of this series, I’ve discussed permutations, generating uniform random points in a triangle, on a sphere, the inversion method, even recycling expensive random bits. In this entry, I will use […]

  10. Darrell says:

    One small optimization to reflecting around x+y=1 to avoid a bit of arithmetic – map the triangle (0,0),(1,0),(1,1) into your target triangle instead and just swap x and y whenever y>x. A bit simpler and no arithmetic at all.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s

%d bloggers like this: