Generate random locations within a triangular domain
Asked Answered
C

3

8

I want to generate x and y having a uniform distribution and limited by [xmin,xmax] and [ymin,ymax]

The points (x,y) should be inside a triangle.

How can I solve such a problem?

Cheerly answered 21/11, 2017 at 9:53 Comment(4)
The requirements here seem underspecified (and possibly contradictory): if x and y are independent and each is uniformly distributed on an interval, they'll cover a rectangle rather than a triangle. How is this triangle specified, and what's its relation to the range [xmin, xmax] and [ymin, ymax]?Haberdashery
Don’t just tell us what you want (especially underspecified). What have you tried? What’s gone wrong? Where are you stuck?Duprey
We probably all guessed that min, max relates to the triangle coordinates. But looking at her other question, I am afraid, we won't get any feedback.Xerosis
See mathworld.wolfram.com/TrianglePointPicking.htmlThroughout
H
15

Here's some code that generates points uniformly on an arbitrary triangle in the plane.

import random
    
def point_on_triangle(pt1, pt2, pt3):
    """
    Random point on the triangle with vertices pt1, pt2 and pt3.
    """
    x, y = sorted([random.random(), random.random()])
    s, t, u = x, y - x, 1 - y
    return (s * pt1[0] + t * pt2[0] + u * pt3[0],
            s * pt1[1] + t * pt2[1] + u * pt3[1])

The idea is to compute a weighted average of the three vertices, with the weights given by a random break of the unit interval [0, 1] into three pieces (uniformly over all such breaks). Here x and y represent the places at which we break the unit interval, and s, t and u are the length of the pieces following that break. We then use s, t and u as the barycentric coordinates of the point in the triangle.

Here's a variant of the above that avoids the need to sort, instead making use of an absolute value call:

def point_on_triangle2(pt1, pt2, pt3):
    """
    Random point on the triangle with vertices pt1, pt2 and pt3.
    """
    x, y = random.random(), random.random()
    q = abs(x - y)
    s, t, u = q, 0.5 * (x + y - q), 1 - 0.5 * (q + x + y)
    return (
        s * pt1[0] + t * pt2[0] + u * pt3[0],
        s * pt1[1] + t * pt2[1] + u * pt3[1],
    )

Here's an example usage that generates 10000 points in a triangle:

pt1 = (1, 1)
pt2 = (2, 4)
pt3 = (5, 2)
points = [point_on_triangle(pt1, pt2, pt3) for _ in range(10000)]

And a plot obtained from the above, demonstrating the uniformity. The plot was generated by this code:

import matplotlib.pyplot as plt
x, y = zip(*points)
plt.scatter(x, y, s=0.1)
plt.show()

Here's the image:

enter image description here

And since you tagged the question with the "numpy" tag, here's a NumPy version that generates multiple samples at once. Note that it uses the matrix multiplication operator @, introduced in Python 3.5 and supported in NumPy >= 1.10. You'll need to replace that with a call to np.dot on older Python or NumPy versions.

import numpy as np

def points_on_triangle(v, n):
    """
    Give n random points uniformly on a triangle.

    The vertices of the triangle are given by the shape
    (2, 3) array *v*: one vertex per row.
    """
    x = np.sort(np.random.rand(2, n), axis=0)
    return np.column_stack([x[0], x[1]-x[0], 1.0-x[1]]) @ v


# Example usage
v = np.array([(1, 1), (2, 4), (5, 2)])
points = points_on_triangle(v, 10000)
Haberdashery answered 21/11, 2017 at 16:57 Comment(2)
Your answer was instrumental in me writing a quick Sierpinski Triangle generator, thought I'd share.Rompers
u could simply be 1-s-t.Jammiejammin
A
9

Ok, time to add another version, I guess. There is known algorithm to sample uniformly in triangle, see the paper by Osada et al. (2002), chapter 4.2 for details.

Python code:

import math
import random

import matplotlib.pyplot as plt

def trisample(A, B, C):
    """
    Given three vertices A, B, C, 
    sample point uniformly in the triangle
    """
    r1 = random.random()
    r2 = random.random()

    s1 = math.sqrt(r1)

    x = A[0] * (1.0 - s1) + B[0] * (1.0 - r2) * s1 + C[0] * r2 * s1
    y = A[1] * (1.0 - s1) + B[1] * (1.0 - r2) * s1 + C[1] * r2 * s1

    return (x, y)

random.seed(312345)
A = (1, 1)
B = (2, 4)
C = (5, 2)
points = [trisample(A, B, C) for _ in range(10000)]

xx, yy = zip(*points)
plt.scatter(xx, yy, s=0.2)
plt.show()

And result looks like

enter image description here

Osada, R., Funkhouser, T., Chazelle, B., & Dobkin, D. (2002). Shape distributions. ACM Transactions on Graphics (TOG), 21(4), 807-832.

Arabel answered 22/11, 2017 at 0:38 Comment(0)
G
2

Uniform on the triangle?

import numpy as np

N = 10 # number of points to create in one go

rvs = np.random.random((N, 2)) # uniform on the unit square
# Now use the fact that the unit square is tiled by the two triangles
# 0 <= y <= x <= 1 and 0 <= x < y <= 1
# which are mapped onto each other (except for the diagonal which has
# probability 0) by swapping x and y.
# We use this map to send all points of the square to the same of the
# two triangles. Because the map preserves areas this will yield 
# uniformly distributed points.
rvs = np.where(rvs[:, 0, None]>rvs[:, 1, None], rvs, rvs[:, ::-1])

Finally, transform the coordinates
xmin, ymin, xmax, ymax = -0.1, 1.1, 2.0, 3.3
rvs = np.array((ymin, xmin)) + rvs*(ymax-ymin, xmax-xmin)

Uniform marginals? The simplest solution would be to uniformly concentrate the mass on the line (ymin, xmin) - (ymax, xmax)

rvs = np.random.random((N,))
rvs = np.c_[ymin + (ymax-ymin)*rvs, xmin + (xmax-xmin)*rvs]

but that is not very interesting, is it?

Garland answered 21/11, 2017 at 10:11 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.