Delaunay Triangulation for Terrain Generation in Unity, Part 2

This post is a companion post to my previous post on using delaunay triangulation in terrain generation. In that post, I covered some of the implementation details involved in taking a Triangle.Net mesh and transforming it into a Unity3D mesh.

In this post, I’ll talk about a couple more things I did with the terrain:

• How to space out the triangulated points better
• How to find the height at an arbitrary point in worldspace

Spacing out the terrain

If you remember, last time we left off with this terrain made from a triangulated set of randomly placed vertices:

If you zoom in a bit, you’ll see that there tons of triangles that are tiny slivers and don’t add much to the overall detail:

It should be pretty obvious why this is happening - randomly spaced points give you no guarantees. What we need is some method to space the points out instead of placing the points randomly.

In amitp’s article on Voronoi terrain generation, he uses Lloyd relaxation to make the polygons more regular. I tried this method, and I couldn’t get it to work efficiently with Triangle.Net - I suspect it would require some extensions to the source code.

Thus, instead of post-processing the points, I chose to create evenly-spaced points upfront. To create evenly spaced points, I used a technique called Poisson Disk Sampling. As in the last article, I won’t go into the full details of how it works (and, in fact, I am still not totally sure how it works). Here’s a decent reference for the implementation.

For my part, I ended up taking this guy’s excellent implementation in c# and adapted it for use in my code. Here’s the result after sampling using poisson disk sampling for sampling points:

This still looks more interesting than a grid-based solution, but is a bit too regular. To remedy the samey-ness of the new terrain, I combined the two sampling methods; I use some points from poisson disk sampling, and some points from random sampling. Here’s the result from that:

You can definitely still see some tiny creases that come from the triangle slivers, but the overall effect is improved from the original.

Querying height of terrain in worldspace

The terrain’s looking pretty nice, but what if we want to put some grass or rocks or trees on it? We’ll need to know the elevation (y) at a given (x, z) coordinate in worldspace. If we had a uniform grid, solving this problem is pretty easy:

1. Divide the worldspace point by the size of the grid to get the square containing the worldspace point.
2. Do some math to figure out the height of the square at the local point.

With a triangulation the steps are similar, but now they become:

1. Find which triangle contains the given point.
2. Do some math to figure out the height of the triangle at the local point.

Some of you may be thinking; why not add a mesh collider to the terrain and use raycasts? To be honest, I can’t come up with anything seriously wrong with this method. I suppose I just find it distasteful to use physics for a problem like this, which is unrelated to physics. So, off we go, with the first step being:

Finding the triangle containing the point

The naive way to solve the problem would be to iterate over all triangles, checking whether the given vertex is contained within the triangle. However, this is monstrously slow.

I went through a couple of thoughts on how to resolve this problem. One way would be to use a classic spatial partitioning structures, such as a quadtree, to bin the triangles. What if we just used a grid, instead of a quadtree?

Well, believe it or not, this problem is somewhat related to the previous problem of generating uniform points. Now that we have evenly spaced points, we should be able to create a grid where, on average, each of the grid spaces contains one vertex.

In fact, if you look closely at the Poission Disk Sampling code, you’ll see that internally, it uses a grid, and does guarantee that each space in that grid contains one point. With this knowledge, we can create a data structure for querying triangles by doing the following:

1. Create a grid to hold the points, the same size as the poisson disk-sampling grid.
2. Bin the vertices into this grid.
3. Create a map of vertices to triangles that they participate in.

Here’s a visual representation I whipped up to (hopefully) clarify the idea:

The green and yellow circles indicate two vertices that would be binned in the (3,2) position in the 2-d array. There is a separate data structure to associate vertices with triangles; the green-shaded triangles are associated with the green vertex, the yellow-shaded triangles with the yellow vertex, and the green-and-yellow-shaded triangles with both. Each vertex is binned in this way.

Even with our randomly-spaced points thrown in, this binning method should be good enough of a spatial partitioning method. To query for triangles, we can do the following:

1. Figure out which bin corresponds to the point we’re looking for.
2. Find all vertices in that bin (and some surrounding ones, for safety).
3. Check which of the participant triangles contain the vertex.

By doing this, we’ve culled the list of possible triangles down from thousands down to tens.

The code for all this is fairly straightforward, but a bit long, so I won’t paste it here. You can find it in the github repo, under the classname TriangleBin.

Obtaining the elevation within the triangle

We now have a triangle, but we still don’t know the actual elevation within that triangle at the given (x, z) point. We can obtain this using a plane-line intersection, where the plane is the one formed by the three points of the triangle (\(t_0\), \(t_1\), and \(t_2\)), and the line is the line containing point \((x, 0, z)\) with direction \((0, 1, 0)\).

For more information on a plane-line intersection, you can refer to the wikipedia article. Here is the set of equations we need:

• \(i\) is the point of intersection.
• \(p_0\) is any point within the plane. We can use any of \(t_0\), \(t_1\), or \(t_2\). I’ll stick with \(t_0\).
• \(n\) is the normal of the plane. It can be calculated using a cross product: \((t_0 - t_1) \times (t_1 - t_2)\).
• \(l_0\) is any point which the line goes through. We can use \((x, 0, z)\).
• \(l\) is the direction of the line, \((0, 1, 0)\).

Note that we only need the y-component of \(i\), and not the other terms, so we can get the elevation \(e = i_y\) by solving the following:

Since \(l_y = 1\) and \(l_{0_y} = 0\), it turns out that \(e = d\). If we simplify the equation for \(d\), we can find the equation in terms of \(x\) and \(z\):

Now that we have both the parts - finding the triangle containing a point, and then finding the elevation within that triangle - we can finally turn all this into some code:

/* Returns the triangle containing the given point. If no triangle was found,
   then null is returned.  The list will contain exactly three point indices. */
public List<int> GetTriangleContainingPoint(Vector2 point) {
    // bin is an instance of TriangleBin
    Triangle triangle = bin.getTriangleForPoint(new Point(point.x, point.y));
    if (triangle == null) {
        return null;
    }

    return new List<int>(new int[] { triangle.vertices[0].id,
                                     triangle.vertices[1].id,
                                     triangle.vertices[2].id });
}

/* Returns a pretty good approximation of the height at a given point in worldspace */
public float GetElevation(float x, float y) {
    // Make sure we're in-bounds
    if (x < 0 || x > xsize || y < 0 || y > ysize) {
        return 0.0f;
    }

    // Get the triangle containing the point
    Vector2 point = new Vector2(x, y);
    List<int> triangle = GetTriangleContainingPoint(point);

    if (triangle == null) {
        // This can happen sometimes because the triangulation does not actually fit entirely
        // within the bounds of the grid - not great error handling on show here, but it works
        return float.MinValue;
    }

    // The implementation of GetPoint3D is covered in the previous post;
    // given a vertex index, it returns the position in local space
    Vector3 p0 = GetPoint3D(triangle[0]);
    Vector3 p1 = GetPoint3D(triangle[1]);
    Vector3 p2 = GetPoint3D(triangle[2]);

    // Here's where we use the line-plane intersection equation derived above
    Vector3 normal = Vector3.Cross(p0 - p1, p1 - p2).normalized;
    float elevation = p0.y + (normal.x * (p0.x - x) + normal.z * (p0.z - y)) / normal.y;

    return elevation;
}

/* Scatters detail meshes randomly within the bounds of the terrain. */
public void ScatterDetailMeshes() {
    for (int i = 0; i < detailMeshesToGenerate; i++)
    {
        // Obtain a random position
        float x = Random.Range(0, xsize);
        float z = Random.Range(0, ysize);
        float elevation = GetElevation(x, z);
        Vector3 position = new Vector3(x, elevation, z);

        if (elevation == float.MinValue) {
            // Value returned when we couldn't find a triangle, just skip this one
            continue;
        }

        // We always want the mesh to remain upright, so only vary the rotation in
        // the x-z plane
        float angle = Random.Range(0, 360.0f);
        Quaternion randomRotation = Quaternion.AngleAxis(angle, Vector3.up);

        Instantiate<Transform>(detailMesh, position, randomRotation, this.transform);
    }
}

Note that there’s one other change we have to make. In the previous post, I scaled the terrain elevation by changing the y scale of the transform. Now that we need the world-space elevation within the code itself, I’ve started scaling the numbers within the elevations array by a passed-in value during the perlin noise sampling.

Now we can call ScatterDetailMeshes within Generate. Make sure to assign detailMesh in the inspector; I’ve assigned a rock mesh to it.

The rocks are a bit hard to see at that distance, but if you come closer, you can see that they are indeed scattered neatly on the surface of the terrain:

That’s it for this tutorial. The github repository containing all the code can be found here. Thanks for reading.

MathJax used for the equations.