Meandering Triangles

An in-depth explanation of the Meandering Triangles algorithm for drawing contour lines.

Estimated reading time: 5 mins

Contour lines, also known as isolines, are lines drawn along along paths of constant elevation. They are often used when drawing maps to indicate the slopes of hills and mountains.

But they can also be used to visualize functions that map (x,y) values to z values, or to find the boundaries of implicit surfaces. One of the common algorithms for finding contour lines is Marching Squares. The algorithm assumes that you have a set of elevation data points spaced out at regular intervals on a grid, and generates contour lines for each square of the grid. A simpler algorithm exists, called Meandering Triangles, and it divides the data points into triangles, rather than squares, to reduce the number of cases the algorithm has to handle. Here, I’ll give a full explanation of how the Meandering Triangles algorithm works, with some examples to show it in action.

Meandering Triangles in Action

First, let’s see how the meandering triangles algorithm looks:

Step 0: Get Some Data

In order to generate contour lines, we need some data whose contours we want to find.

import math
# some arbitrary function that maps (x,y) to an elevation between 0 and 1
def elevation_function(x,y):
    return 1/(2+math.sin(2*math.sqrt(x*x + y*y))) * (.75+.5*math.sin(x*2))

elevation_data = dict()
WIDTH, HEIGHT = 100, 100
SPACING = 1
for x in range(0,width, SPACING):
    for y in range(0,height, SPACING):
        elevation_data[(x,y)] = elevation_function(x,y)

Step 1: Triangulate

The first step is to divide the regularly spaced grid of elevation data into triangles. In this example, I divide each set of 4 adjacent points into the top-left and bottom-right halves of a square, which forms two triangles.

import collections
Triangle = collections.namedtuple("Triangle", "v1 v2 v3")
triangles = []
for x in range(0,width-1, SPACING):
    for y in range(0,height-1, SPACING):
        t1 = Triangle((x,y), (x+SPACING,y), (x,y+SPACING))
        triangles.append(t1)
        t2 = Triangle((x+SPACING,y), (x,y+SPACING), (x+SPACING,y+SPACING))
        triangles.append(t2)

Step 2: Find Triangles That Intersect Contour Lines

If we assume the elevation data is sampled at a high enough resolution, then each triangle is pretty close to a small, nearly-flat bit of terrain. If we’re drawing a contour line at elevation = 0.5, then the triangle will have a contour line going through it if some part of that triangle is below 0.5, and some part of it is above 0.5. The same is true for any threshold, but let’s assume we want to draw the contour line for 0.5 in this example. There are four possible cases here:

• All the vertices are below the threshold (no contour line)
• All the vertices are above the threshold (no contour line)
• Two of the vertices are below, and one is above the threshold (has a contour line)
• One of the vertices is below, and two are above the threshold (has a contour line)

Here’s the code:

threshold = 0.5
Edge = collections.namedtuple("Edge", "e1 e2")
contour_segments = []
for triangle in triangles:
    below = [v for v in triangle if elevation_data[v] < threshold]
    above = [v for v in triangle if elevation_data[v] >= threshold]
    # All above or all below means no contour line here
    if len(below) == 0 or len(above) == 0:
        continue
    # We have a contour line, let's find it

Step 3: Find the Contour Line Segments

Once we know that a contour line passes through a triangle, we also know that it passes between the vertices that are above the threshold and those that are below the threshold. Since it’s a triangle, we have 1 vertex on one side, and 2 on the other side, and the contour line will pass through 2 of the 3 edges of the triangle. Where along that edge the contour line will cross can be determined by linearly interpolating between the triangle vertices to find where the threshold’s exact value should fall.

threshold = 0.5
Edge = collections.namedtuple("Edge", "e1 e2")
contour_segments = []
for triangle in triangles:
    below = [v for v in triangle if elevation_data[v] < threshold]
    above = [v for v in triangle if elevation_data[v] >= threshold]
    # All above or all below means no contour line here
    if len(below) == 0 or len(above) == 0:
        continue
    # We have a contour line, let's find it
    minority = above if len(above) < len(below) else below
    majority = above if len(above) > len(below) else below

    contour_points = []
    crossed_edges = (Edge(minority[0], majority[0]),
                     Edge(minority[0], majority[1]))
    for triangle_edge in crossed_edges:
        # how_far is a number between 0 and 1 indicating what percent
        # of the way along the edge (e1,e2) the crossing point is
        e1, e2 = triangle_edge.e1, triangle_edge.e2
        how_far = ((threshold - elevation_data[e2])
                   / (elevation_data[e1] - elevation_data[e2]))
        crossing_point = (
            how_far * e1[0] + (1-how_far) * e2[0],
            how_far * e1[1] + (1-how_far) * e2[1])
        contour_points.append(crossing_point)

    contour_segments.append(Edge(contour_points[0], contour_points[1]))

Step 4: Joining Up

All that’s left to do is join up all the adjacent line segments into lines. My implementation builds up contour lines by adding segments to the head and tail of a line until it runs out of connected segments.

unused_segments = set(contour_segments)
segments_by_point = collections.defaultdict(set)
for segment in contour_segments:
    segments_by_point[segment.e1].add(segment)
    segments_by_point[segment.e2].add(segment)

contour_lines = []
while unused_segments:
    # Start with a random segment
    line = collections.deque(unused_segments.pop())
    while True:
        tail_candidates = segments_by_point[line[-1]] & unused_segments
        if tail_candidates:
            tail = tail_candidates.pop()
            line.append(tail.e1 if tail.e2 == line[-1] else tail.e2)
            unused_segments.remove(tail)
        head_candidates = segments_by_point[line[0]] & unused_segments
        if head_candidates:
            head = head_candidates.pop()
            line.appendleft(head.e1 if head.e2 == line[0] else head.e2)
            unused_segments.remove(head)
        if not tail_candidates and not head_candidates:
            # There are no more segments touching this line,
            # so we're done with it.
            contour_lines.append(list(line))
            break

And that’s it! The example code I’ve shown here assumes a single fixed threshold, but it’s easy to run multiple passes to find contour lines at multiple different thresholds. I think this algorithm is quite elegant in its simplicity, and it runs quickly on arbitrary data. There are certainly better algorithms if you care about the level of detail and have an elevation function that is slow to compute, but this algorithm makes an excellent first draft for finding contour lines.

You can get the full code here.