A Brief Guide to the R “Triangulation” Package

Introduction

A triangle mesh is a type of polygon mesh constructed over some 2D spaces. It comprises a set of triangles that are connected by their common edges or vertices covering the space of interest. Triangulation is widely used in various areas, such as Computer Graphics, Geography, Environmental Science, Medical Studies, and Statistics. There are many useful toolboxes that can be used for generating a 2D triangle mesh, such as (i) the “simple mesh generator” based on MATLAB DistMesh, (ii) the “two-dimensional quality mesh generator and Delaunay triangulator” based on C language Triangle, and (iii) the R package Triangulation. This blog will give you a brief introduction to the Triangulation R package.

The Triangulation R package is a handy tool for generating and visualizing triangle meshes of a 2D domain of arbitrary shape. This package contains two main functions, (i) the triangle mesh generator, TriMesh(), and (ii) the visualization function for triangulations, TriPlot(). It also includes the datasets which are necessary for multiple toy examples, such as the boundary of a horseshoe domain (horseshoe), the exterior boundary of the mainland of the United States (USbb), and boundaries of several irregular shapes (shape and weird). For detailed information on each of the data (BMP, horseshoe, shape, weird, USbb) and function (TriMesh() and TriPlot()) available in Triangulation, please refer to the R documentation.

There are several Delaunay triangulation methods, including (i) exact Delaunay triangulation, (ii) constrained Delaunay triangulation, (iii) conforming Delaunay triangulation, and (iv) constrained conforming Delaunay triangulation. The exact Delaunay triangulation is constructed based on the position of a given set of vertices while not regarding much on the ways of their connection. It tries to maximize the minimum angle within the triangulation. Unlike the exact Delaunay triangulation, the constrained Delaunay triangulation could force some specifically required segments into the triangulation as edges. Besides the given set of vertices, a conforming Delaunay triangulation usually inserts additional vertices when forming a Delaunay triangulation to allow the existence of the segments in the mesh, meeting constraints on the minimum angle and maximum triangle area while maintaining the Delaunay property. The function TriMesh() in the Triangulation package focuses on the conforming Delaunay triangulation to triangulate the domain.

How to install

The Triangulation package is currently available on GitHub. To install the package directly from GitHub, first we need another R package, devtools. It can be installed as follows:

Now, the Triangulation package can be installed using the following code:

Next, we need to load the Triangulation package to use:

How to use

This section will discuss the practical implementation of the Triangulation package on a few datasets that are supposed to cover different types of domains that we usually come across in real-life scenarios.

Example for Rectangular Domain

First, we consider a basic rectangular domain(without any hole) defined on \([0, 2] \times [0, 1]\)

The TriMesh() function can then be used to triangulate this rectangular domain. We have the following three input arguments when use the TriMesh() function.

  • Pt: An \(N\) (number of bountary vertices) \(\times 2\) matrix indicates the outer boundary points of the 2D domain.
  • n: An integer parameter controlling the fineness of the triangulation.
  • H: A list of vertices denotes the inner boundary points (when holes are present), and the default value is NULL.

## The number of vertices in this triangulation is  231  and the number of triangles is 400

The output of TriMesh() contain two components:

  • V: An \(N_V\) (number of vertices) \(\times 2\) matrix, with \(i\)th row being the Cartesian coordinates for the \(i\)th vertex in the triangulation.

  • Tr: A \(K\) (number of triangles) \(\times 3\) matrix, with each row representing one triangle. All the elements in a row, are the integers that stand for the indices of vertices in \(V\).

The argument n in the TriMesh() function can be used to fine-tune the mesh construction. A larger value of n usually results in a finer triangulation. For example, if we increase n to 20 in the previous example, we will have the following finer triangulation of the same rectangular domain.

## The number of vertices in this new triangulation is  861  and the number of triangles is 1600

With finer triangulation (i.e., with higher n), the number of vertices and the number of triangles generated also increases.

Even though we can automatically get the triangle mesh plot from TriMesh(), we can also obtain the same plot using the following commands.

Examples for Complex Domains

1. Horseshoe domain

The horseshoe domain reflects an example that the domain is almost divided into two parts by a strict concavity with a very narrow gap between the two segments of the domain.

Specifically, the boundary of the domain is shown in the following figure.

And the first few boundary points (coordinates) of the Horseshoe domain are

##              V1           V2
## 1 -9.000000e-01 1.102146e-16
## 2 -7.534498e-01 4.922533e-01
## 3 -4.283527e-01 7.915264e-01
## 4  5.510729e-17 9.000000e-01
## 5  3.000000e+00 9.000000e-01
## 6  3.240297e+00 8.197771e-01

Using the TriMesh() function, we can obtain a triangulation of the horse domain as follows.

2. Irregular domain with holes:

The BMP data contains a list of coordinates of the boundary vertices of Bandiagara, a small town in Mali. It represents a case in our real life in which the closed bounded domain contains some holes.

By using the H argument in our TriMesh() function, we can construct a triangulation of the domain:

3. US map

Now let’s try another real example and build a triangle mesh for the mainland of the United States. The US boundary can be obtained as follows:

##      x.coord  y.coord
## 1 -122.59593 49.05139
## 2 -114.97487 49.05139
## 3 -106.02840 49.05139
## 4  -98.87123 49.05139
## 5  -94.89502 49.44901
## 6  -90.72001 48.18988

Again, using TriMesh(), we can get the triangulation for US:

Few notes to consider

  • The parameter n in the TriMesh() function should be chosen in such a way that
    • the triangles that are generated are more or less uniform in size;
    • the angles of the triangles are not too acute.
  • To apply TriMesh(), we need to have the polygonal vertices that represent the outer boundary of the domain under consideration and the boundary of the holes within the domain.