Uniform sampling of a triangle by dividing it into smaller parts?

user_2942145 0
Viewed 1630

To sample a triangle ABC uniformly, I can use the following formula:

P = (1 - sqrt(r1)) * A + (sqrt(r1)*(1 - r2)) * B + (r2*sqrt(r1)) * C

where r1 and r2 are random numbers between 0 and 1. The more samples you take, the better. But what if I want to get a better distribution, while keeping then number of samples low?

For example if I had a square, I can implicitly divide it into an N x N grid and generate a random sample inside the smaller grid squares. Like this:

float u = (x + rnd(seed)) / width;
float v = (y + rnd(seed)) / height;

The point is I force the sampling to cover the entire grid at a lower sample resolution.

How can I achieve this with a triangle? The only way I can think of is to explicitly subdivide it into a number of triangles using a library like Triangle. But is there a way to do this implicitly like with a square, without having to actually divide the triangle?

triangulation random
2 Answers

OK, I had some thoughts and believe using quasirandom numbers could improve "uniformity" of the points-in-the-triangle coverage without doing subdivision into smaller triangles. Quasirandom sampling using Sobol sequences could provide a lot better coverage as seen in the Wiki article.

Here is 200 points in triangle using standard RNG (whatever it is in Python) enter image description here

And here is picture with 200 points sampled from Sobol 2D sequence

enter image description here

Looks a lot better to me. Python code to play with

import os
import math
import random

import numpy as np
import matplotlib.pyplot as plt

import sobol_seq

def trisample(A, B, C, r1, r2):
    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)

if __name__ == "__main__":
    N = 200

    A = (0.0, 0.0)
    B = (1.0, 0.0)
    C = (0.5, 1.0)

    seed = 1
    xx = list()
    yy = list()
    random.seed(312345)
    for k in range(0, N):

        pts, seed = sobol_seq.i4_sobol(2, seed)
        r1 = pts[0]
        r2 = pts[1]

        # uncomment if you want standard rng
        #r1 = random.random()
        #r2 = random.random()

        pt = trisample(A, B, C, r1, r2)
        xx.append(pt[0])
        yy.append(pt[1])

    plt.scatter(xx,  yy)
    plt.show()

    sys.exit(0)
user_40446960
Related