Triangulation with Python
Ever wondered how to connect multiple 3-dimensional points into a surface? The answer is triangulation. We’re going to look at what it is and how to use it in python.
Before we start
The following context will be used throughout this post.
The code from the post is also on this repository.
What is triangulation?
Say we have a set of points in 3d space, which describe a surface. But for reasons regarding optimization or lack of data, we can’t expand those points to form a continuous surface. Enter triangulation. What if we describe a set of triangles with our points as vertices and draw them instead? That sounds perfect!
But how do we define those triangles? Luckily, there exists an algorithm called Delaunay triangulation, that can define them for us. Though it does that only for sets of 2d points, as in 3d there’s no one way of triangulating a given set of points (different triangulations may produce different surfaces).
The algorithm will make each triangle in a way that its vertices are as close to one another as possible. The reason is, that we can’t have a triangle drawn over three distant points that may have other points in-between them. Triangles must all be glued up together side-by-side, without any one of them containing another.
Using the algorithm, it is easy to triangulate 3d surfaces that are described as points (x, y, f(x,y)), where f is some function. That means every pair (x, y) has exactly one point z (height) assigned to it. Our triangulation of the points (x, y, f(x,y)) is defined via the triangulation of projections (x, y), as they completely describe the closeness of points on our surface.
A slight problem arises when we have a more complex 3d surface, like a sphere. Then we can’t project the points onto a plane, as that no longer preserves the information about what points were close to one another. For example points (0,0,1) and (0,0,-1) lie on opposite sides of the unit sphere, but would both project to the same point (0,0) on the plane. In that case, we have to map to our surface from a triangulated plane. This is the definition of parametrizing a surface.
Note that we also did this before, as we mapped (x,y) to (x, y, f(x,y)).
See an example of how we fold a piece of paper into a torus. This folding can be viewed as a parametrization. It should be clear, how we get a triangulation of a torus from the triangulation of points on the paper (we leave the triangles alone during the folding process).
Triangulating a face
Here we’ll triangulate this set of 3d points that describe a face. In our case, the face is described by points of the form (x, y, f(x,y)), so we can use the easy approach of triangulating their projections.
We also added cm.inferno color map to visualize the height of points. After running the code, we get something like this.
Triangulating an ellipsoid
This example is inspired by my statistics homework. We were given the task of linearly regressing some data and calculating the confidence ellipsoid of regression parameters. That is, we had to draw an ellipsoid, for which the probability of those parameters being inside, had been large.
The moral of the story is that drawing ellipsoids is pretty useful. And you’d think there’d be some matplotlib function, allowing you to do this. Well, you’d be wrong. I had to start from scratch and make use of triangulation.
A sphere
To not bore you with exact details of how ellipsoids are drawn right at the beginning, let’s first learn how to draw a (unit) sphere. We’ll see that from then on drawing an ellipsoid is pretty easy.
We need to parametrize our sphere. The most natural parametrization is
(cos(u) sin(v), sin(u) sin(v), cos(v)),
where u goes from 0 to 2π and v goes from 0 to π.
If you’re interested, this is derived from the fact that given any two normed and orthogonal vectors q and p, the parametrization cos(t) q + sin(t) p (where t goes from 0 to 2π) represents a unit circle, spanned over q and p. To draw a sphere, we draw (half) a circle through each pair of vectors (0, 0, 1) and (cos(u), sin(u), 0), where u goes from 0 to 2π.
Let’s implement this.
Increasing the parameter k defines more points, thus smaller triangles and thus a smoother sphere. Running the code produces something like this:
You can also draw parts of the sphere and observe its interior, by limiting the domain to e.g. [0, 2𝝅] x [0, 𝝅/2].
An ellipsoid
All that we have to do to get an ellipsoid is pick a 3x3 matrix and apply it to each point on the sphere. The matrix defines the axes of an ellipsoid and the stretching along those axes. In the context of statistics, this is closely related to the covariance matrix.
The matrix we used was the one from my homework.
Triangulating two other surfaces
Remarks
There are algorithms besides Delaunay triangulation. Some can triangulate actual 3d points, so we don’t need to bother with parametrizations. For example, ConvexHull can handle this for convex objects.
