Triangular Grid Sampling

Last summer, I worked out a really cool algorithm for sampling points on a triangular grid. The core idea is shockingly elegant, and I've always wanted to explain it in detail.

In my use case, I wanted to numerically approximate an integral over the 2D area of a triangle. The simplest way is to use a Riemann sum, and I needed a grid of points to evaluate my function over.

\int \mathtt{my\_ func(r)}\; dA \approx \sum_i \mathtt{my\_ func(r_i)}\; \delta A

I found a beautifully simple approach to creating a triangular grid that relies on the self-tiling nature of triangles. You can split any triangle into four smaller, identical triangles by connecting the midpoints of its sides. Three of these triangles have the same orientation as the parent, while the central one is flipped upside down.

This self-similarity is the key. If we can generate a grid for the large triangle, we can simply scale and translate that same grid pattern to fit into each of the four sub-triangles.

We can fill the three corners of the large triangle with the shrunken down smaller triangle, but we need to flip the smaller triangle in order to fit it in the middle slot.

You can imagine doing this same process again. We scale down our samples...

and then translate, copy and flip.

Isn't this so cool? The final code is delightfully simple.

fn samples(corners: [Position; 3], iters: usize) -> Vec<Position> {
    let centroid: Position = corners.iter().sum() / 3.0;
    let mut samples = vec![centroid];

    for _ in 0..iters {
        let mut new_samples = Vec::with_capacity(samples.len() * 4);

        // scale all points down by a factor of 1/2
        // creates the 'shrunken down' green dots as seen earlier!
        samples.iter_mut().for_each(|x| *x /= 2.0);

        // fill in the three corners
        for c in corners {
            new_samples.extend(samples.iter()
                .map(|half_sample| half_sample + c / 2.0))
        }

        // fill in center triangle
        new_samples.extend(
            samples
                .iter()
                .map(|half_sample| 1.5 * centroid - half_sample),
        );
        samples = new_samples;
    }

    samples
}

This algorithm gaurantees that all points are inside the triangle by maintaining the convex hull property. We can prove this with barycentric coordinates. Any point \mathbf r inside the triangle can be written as \mathbf r= a \mathbf c_1 + b \mathbf c_2 + c \mathbf c_3 where a+b+c = 1 and a,b,c \in [1,0]. In both the corner sub-triangle (half_sample + c / 2.0) and in the central sub-triangle (1.5 * centroid - half_sample), the resulting barycentric coordinates are also non-negative and sum to one.

From a performance perspective, this code is excellent! The main loop is branch free and SIMD friendly, generally only needing addition and subtraction. We linearly iterate over samples and linearly write to new_samples in a very cache friendly manner.