** Most two dimensional quasirandom methods focus on sampling over a unit square. However, quadrature over the triangle is also very important in computer graphics. Therefore, I describe a direct construction method for a point set or sequence that evenly covers an arbitrary shaped triangle. **** ****I first show how to create an evenly distributed (finite) point set over the triangle, and then present a method of generalising it to an open (infinite) quasirandom low discrepancy point sequence over the triangle. Finally, I briefly discuss how this technique can be applied to other geometries.**

**Outline**

Low discrepancy sequences that evenly sample/cover a square have been extensively studied for nearly a hundred years. Most of these quasirandom sequences can be extended to rectangles by simple dilation without much deterioration of discrepancy.

However, this post looks at the interesting and important extension of low discrepancy sequences over an arbitrary triangle.

As far as I can deduce, very little work has been done on constructing evenly distributed points sets and sequences over triangles. Notable efforts in recent years by Basu [2014] and Pillands, [2005] and Brandolini [2013] represent the main papers – if not the only papers on this subject..

These authors generally take a very theoretical and analytical approach to this issue and almost solely focus on the unit regular triangle. In contrast, this post is mainly intended to simply aid in the development of some practical techniques for graphics rendering.

This post describes a simple and direct method that requires no rejection, dropping or recursion; and most importantly it can be applied to triangles of arbitrary shape and size.

This post has four sections

- The basic Fibonacci lattice
- Evenly distributing a set of $N$ points in a triangle.
- Constructing an open (infinite) set of evenly distributed points in a triangle.
- Further generalizations

**1. The Fibonacci Lattice**

You might think it is easy to place 100 points in a square such that the minimum distance between neighboring points is as large possible. But what if you had to place 13 points? 47? … or how about 2019 points?!

I doubt that you could do better than those shown in figure 2. But don’t worry. The not-so-hidden secret to constructing these three point distributions (and any others you may need) is that there is an amazingly simple algorithm that can produce them all.

**Enter the Fibonacci lattice.**

We define the Fibonacci lattice as of $N$ points where: $ \pmb{p}_i(x,y) = \left( r_1(i,r_2(i) \right) $ where

$$ r_1(i) = \left( \frac{i-0.5}{N} \right) ;\quad \textrm{and } r_2(i) = \left\{ \frac{i}{\phi} \right\} \quad \textrm{for } \; i=1,2,3,…, N. \tag{1}$$

and $\phi =(1+\sqrt{5})/2 \simeq 1.618033988$ is the golden ratio. Note that the expression $\{ x \} $ denotes the fractional part of $x$. In computing, this fractional operator is usually written and implemented as $ x \; \%1$ (pronounced ‘mod 1’), and as this post as mainly geared towards computing applications we shall use this notation from here onwards. The Fibonacci lattice has many beautiful properties.

This square lattice can be transformed into a disk/spiral with evenly-spaced points, via $ (r,\theta) \leftarrow (\sqrt{r_1}, 2\pi r_1)$, famously known as the Fibonacci spiral. It is very frequently cited that this is the structure of many flowers such as the sunflower due to phyllotaxis.

Also, it can be used as the basis for evenly mapping a set of points on to the surface of the sphere. For an interactive demo of this, see here. (And some small ways that this algorithm be improved are discussed here.)

Also because this is topologically equivalent to a torus (ie a donut), you can tessellate this square tile and it will be seamless tile.

However, one of the key properties that we will use, is that the Fibonacci lattice is what some call ‘scale invariance’. Notice in figure 4 that as the original square lattice is dilated, the evenness remains virtually perfect!

Thus, for any $N$, the point construction: $ \pmb{p}_i(x,y) = (r_1w , r_2h ) $ results in an evenly distributing point set over a rectangle of dimension $(w,h)$.

# 2. Evenly Distributing a set of points in an arbitrary triangle

### Right-angled Triangles

To restrict the placement of points to a region in the shape of a right-angled triangle, we use the following transformation. $$\pmb{p}_i(x,y) = (\sqrt{r_1}, r_2 \sqrt{r_1} ) \tag{2}$$

For very similar reasons to the Fibonacci spiral transformation, the inclusion of the square root is because at any point, the vertical height of the triangle is proportional to the horizontal distance traversed, and since we need the transformation to be area-preserving, we require taking the square root.

And again, this can be dilated in the horizontal or vertical directions. But this time the evenness varies, especially when dilated vertically.

### Arbitrary Triangles

To progress from a right-angled triangle to any shaped triangle now requires introducing non-orthogonal vectors (previously $\pmb{AC}$ was orthogonal to $\pmb{BC}$).

$$\pmb{p} = \pmb{A} + \sqrt{r_1} (\pmb{C} – \pmb{A}) + r_2 \sqrt{r_1} (\pmb{B}-\pmb{C}) \tag{3}$$

where $\pmb{A}, \pmb{B},\pmb{C}$ are the 2-dimensional vectors (coordinates) of the vertices.

That is, equation (3) says that we start at $A$, and move a fraction related to $r_1$ along the side $AC$, and then a fraction related to $r_2$ in the direction parallel to $\pmb{CB}$.

It is not assured that applying such an oblique transformation will preserve discrepancy!

However, fortunately, the figure below shows the result of applying equation (2) to equation (3), for six random triangles and varying $N$, generally preserves the low discrepancy properties very well!

**Side note:** there are two other well known methods to map the two variables $r_1,r_2$ to the inside a triangle: (i) the Parallelogram method and (ii) the Kraemer Method. The acceptance/rejection variation of the Parallelogram method is obviously not applicable here, as we need to keep the point count at exactly $N$. The second variation which involves reflection is equivalent to the Kraemer method. Although these preserve point count, and so could be used, empirical tests show that the reflection/sorting causes a very notable reduction in the evenness / discrepancy of the point set.]

This oblique transform certainly introduces more aliasing / banding than simple horizontal dilations. Techniques to reduce the aliasing will be discussed in a bit more in the next section.

One potentially limiting factor of using the Fibonacci lattice, is that the number of points, $N$ must be known in advance. If after placing $N$ points, you decide you need to place an additional point, you will need to delete all the previous points, and calculate all $N+1$ points from scratch.

Another potential limitation is that the points are placed in the triangle in scan-line order, starting at $\pmb{A}$ and moving towards the $\pmb{BC}$, with the front always parallel to $\pmb{BC}$.

The next section shows how we can resolve both these limitations.

# 3. Evenly Distributing an open sequence of points in an arbitrary triangle.

### Quasirandom Sequences

To evenly distribute an open (infinite) sequence of points inside an arbitrary triangle requires equation (1) to have no explicit dependence on $N$.

One of the simplest ways of constructing an evenly distributing an open sequence of points inside a triangle is to replace the Fibonacci lattice (as defined in equation 1), with a two dimensional low discrepancy quasirandom sequence, and then again apply it to equation 2 or 3.

In this regard, there are several options available to us, including the simple Halton sequence to the more sophisticated Sobol’ sequence. It turns out that depending on how you define ‘evenness’, each of these particular sequences have their advantages and disadvantages. For example, Halton sequences tend to be more useful when placing objects in a region, as the local distance metrics such as nearest neighbor are well optimized. In contrast, the Sobol’ sequence tends to have more ‘clumping’, however, its global point distribution is very even, and so it has excellent quadrature properties.

There is another sequence I love to use, which has both excellent local as well as global properties. This is the $R_2$ sequence, as described in detail in my other post, “The Unreasonable effectiveness of Quasirandom Sequences”.

In brief, we define the infinite two-dimensional sequence $R_2$, such that

$$ \pmb{t}_n = \left\{ n \pmb{\alpha} \right\} \quad n=1,2,3,… $$

$$ \textrm{where} \quad \pmb{\alpha} =\left(\frac{1}{g}, \frac{1}{g^2} \right), $$

$$ \textrm{and} \; g \simeq 1.32471795572 \tag{4}$$

More details on this special value of $g$, which is often called the ‘plastic constant’ can be found on Wikipedia or Mathworld. The first n=150 terms of this $R_2$ sequence are shown in figure 8. Notice that this sequence (like almost all standard low discrepancy sequences) progressively fills the square.

### Anti-aliasing

As mentioned before, one can see in figures 6 and 7 that in certain situations aliasing/banding is clearly visible. This is unfortunately also present (but to a slightly lesser extent) when applying the $R_2$ sequence to the triangular-transform (see figure 1).

**Improvement #1**.

It frequently occurs with triangles that have small angles, and/or triangles that are almost isosceles. It can be reduced by noting that changing the order of the points $A,B,C$ will change the results of equation (2).

Thus with the use of a metric that quantifies or approximates, the aliasing, one could evaluate all six different ways to order the the vertices $ABC, ACB, BAC, BCA, CAB, CBA$, and then from this information, select the preferred version. I have found that the appearance of aliasing is highly dependent on, not only the geometry of the triangle, but also on the particular number of points being placed.

Thus, without testing all six combinations each time, it is a highly non-trivial task to perfectly optimize this selection process. However, I have empirically found that if the vertices are selected in ascending order of the size of vertex angles consistently gives close to minimal aliasing. That is, the recommended implementation applies equation (2) based on the protocol that $A$, $B$ and $C$ are selected in ascending order of angles. (So, $A$ is the smallest angle, and $C$ is the largest angle). All examples and illustrations in this post adopt this heuristic.

**Improvement #2.**

Also one will notice that if instead of the $R_2$ sequence, you use a different low discrepancy sequence such as Halton or Sobol, you will notice that the most regular point distributions are produced using the Fibonacci Lattice, and then by the R_2 an then Halton and lastly then Sobol. However, often the side-effect of such regularity is that these point sequences that are most susceptible to aliasing. That is, aliasing is most visible in the Fibonacci sequence, followed by the $R_2$ sequence, with the Halton and Sobol sequences rarely showing aliasing.

Thus, you need to find a balance between these two tensions of higher regularity but low aliasing.

The solution that I have found works the best is to use a jittered version of $R_2$. For more details on this, see “A simple method to create isotropic quasirandom blue noise“. In this instance, a jitter parameter value of $\lambda \simeq 0.5$ seems optimal. With this value, the regularity Its regularity is almost as good as the Fibonacci sequence, and yet it shows very minimal aliasing.

Thus, combining a jittered version of the $R_2$ sequence with Equation 3, produces the point distributions as shown in figure 8.

Similarly figure 1 also shows this sequence, and demonstrates how it progressively and evenly covers arbitrary triangles.

(Note that to aid comparison, the shapes in figure 9 are the same shaped triangles that are used in figure 7.)

*Figure 9. illustration of the quasirandom low discrepancy sequence mapped to arbitrary triangles and for varying values of $N$. See also figure 1.*

# 4. Further Generalizations

### Higher dimensions

This post has focused on two dimensions, however, equation (2) is amenable to generalisation for higher dimensions via formula such as Grimme, and Smith. For equation (3), the $R_2$ sequence is also fully and naturally generalisable to higher dimensions, as described in detail at “The Unreasonable Effectiveness of Quasirandom Sequences”. Therefore, this method can almost certainly be applied to evenly distributing points over any arbitrary sized and shaped $d-$simplex.

For a generalized Fibonacci lattice, select the first coordinate to vary from 0 to 1, and the other $d-1$ variates to be based on the $d-1$ dimensional quasirandom sequence. In contrast, for the generalized open sequence, simply use the expression for the $d-$dimensional quasirandom sequence.

### Other Geometric Shapes

In theory, this technique can be used easily adapted for other shapes, too.

For example, equation (2) can be generalized to other shapes where the basis vectors are at right-angles.

Briefly, for $x \in [0,1]$, we define $f(x) = 2x$, which is the height of the unit-area triangle.

We then calculate the cumulative area of the region from the origin up to this point. That is, the integral of

$$g(x) = \int f(x) dx = x^2$$

We then need to calculate the inverse cumulative distribution (area).

That is $h(x) = \sqrt{x}$.

Combining all of these, we get the more formal version of equation (1).

$$\pmb{p}_i(x,y) = \left( h(r_1), r_2 f ( h(r_1)) \right) \tag{5}$$

With this equation in hand, we can replace $f(x)=x$, with almost any function we wish, provided we can calculate the inverse cumulative distribution function. In a similar vein to aliasing that appears in figure 6, if you select a function with high peaks, it may result in noticeable aliasing.

A example, I now give is the low discrepancy sequence $R_2$ mapped onto the the classic Gaussian function, is shown in figure 9.

Even further effects can be obtained by leveraging the power of equation (3) which mimics the effect of shearing. A shear is equivalent to thinking of a the graph of a function, where the two axes are not orthogonal.

My name is Martin Roberts. I have a PhD in theoretical physics. I love maths and computing. I’m open to new opportunities – consulting, contract or full-time – so let’s have a chat on how we can work together! Come follow me on Twitter: @Techsparx! My other contact details can be foundhere.