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 slope 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.

First, let’s see how it 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)

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.

## Full Code

''' Copyright © 2017 Bruce Hill <bruce@bruce-hill.com> This work is free. You can redistribute it and/or modify it under the terms of the Do What The Fuck You Want To Public License, Version 2, as published by Sam Hocevar. See http://www.wtfpl.net/ for more details. ''' import math import collections Triangle = collections.namedtuple("Triangle", "v1 v2 v3") Edge = collections.namedtuple("Edge", "e1 e2") def find_contours(elevation_function, xmin, xmax, ymin, ymax, spacing=1, elevation=0.5): ''' Return a list of contour lines that have a constant elevation near the specified value. Lines will be truncated to the sepcified x/y range, and the elevation function will be sampled at even intervals determined by the spacing parameter. ''' elevation_data = dict() for x in range(xmin, xmax+1, spacing): for y in range(ymin, ymax+1, spacing): elevation_data[(x,y)] = elevation_function(x,y) triangles = [] for x in range(xmin, xmax, spacing): for y in range(ymin, ymax, spacing): triangles.append(Triangle((x,y), (x+spacing,y), (x,y+spacing))) triangles.append(Triangle((x+spacing,y), (x,y+spacing), (x+spacing,y+spacing))) contour_segments = [] for triangle in triangles: below = [v for v in triangle if elevation_data[v] < elevation] above = [v for v in triangle if elevation_data[v] >= elevation] # 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 how_far = ((elevation - elevation_data[triangle_edge.e2]) / (elevation_data[triangle_edge.e1] - elevation_data[triangle_edge.e2])) crossing_point = ( how_far * triangle_edge.e1[0] + (1-how_far) * triangle_edge.e2[0], how_far * triangle_edge.e1[1] + (1-how_far) * triangle_edge.e2[1]) contour_points.append(crossing_point) contour_segments.append(Edge(contour_points[0], contour_points[1])) 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]].intersection(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]].intersection(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 return contour_lines if __name__ == "__main__": # Run a test case with 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)) print(find_contours(elevation_function, 0,100, 0,100, spacing=1, elevation=0.5))