Generating Random Points in a Tetrahedron  

Introduction

The problem of picking a random point into a tetrahedron with uniform distribution can arise in various situation, typically in Scientific Visualization when managing unstructured datasets represented by tetrahedral meshes. In some applications there can be the need of creating more than 10000000 samples, so an efficient technique to generate a random sample inside a tetrahedron can be very useful.

Sampling a tetrahedron

The method proposed here is a generalization of one of the techniques used by Turk to generate a random point inside a triangle. The main idea of the 2D technique is to generate a random point in a parallelogram and reflecting it around the center of the parallelogram.

Let s, t and u be three numbers chosen form a uniform distribution of random numbers in the interval [0,1] and A,B and C three 3D vectors; the point sA+tB+uC identifies a random point with uniform distribution inside the 3D skewed parallelepiped defined by A,B and C. To simplify the discussion we will focus only onto s, t and u and how how to dissect and fold the cubic parametric space defined by s, t and u into a tetrahedron. It worth to be noted that while the area of a triangle is 1/2 of the area of the corresponding parallelogram, the volume of a tetrahedron is 1/6 of the volume of the corresponding parallelepiped, therefore we will dissect and fold it in two and the result again in three parts in order to obtain the 2*3=6 tetrahedra. The first step is to cut the cube with the plane s+t=1 into two triangular prisms of equal volume, and then folding all the points falling beyond the plane s+t=1 (i.e. in the upper prism) into the lower prism as shown in figure:

  this can be done by calculating the new (s,t,u) values as shown in the equation:

The second step is to cut and fold the resulting triangular prism with the two planes t+u=1 and s+t+u=1. This dissection identifies the three equal volume tetrahedra shown in figure

 the folding of the triangular prism into the first tetrahedron can be done by calculating the new (s,t,u) values as shown in equation:

These two steps can be straightforwardly implemented in the C++ function that takes in input four 3D points describing the vertexes of a tetrahedron and returns a random point with uniform distribution inside it. We assume that the Point3 class has operators overloaded with the obvious meaning of product between vectors and scalar, sum between vectors, etc.

Point3 Pick(Point3 v0, Point3 v1, Point3 v2, Point3 v3) {

double s = drand();
double t = drand();
double u = drand();
if(s+t>1.0) { // cut'n fold the cube into a prism

s = 1.0 - s;
t = 1.0 - t;

}
if(t+u>1.0) { // cut'n fold the prism into a tetrahedron

double tmp = u;
u = 1.0 - s - t;
t = 1.0 - tmp;

} else if(s+t+u>1.0) {

double tmp = u;
u = s + t + u - 1.0;
s = 1 - t - tmp;

}
a=1-s-t-u; // a,s,t,u are the barycentric coordinates of the random point.
return v0*a + v1*s + v2*t + v3*u;

}