Every 2x2 block of pixels in the binary image forms a contouring cell, so the whole image is represented by a grid of such cells. Note that this contouring grid is one cell smaller in each direction than the original 2D field. For each cell in the contouring grid:
The contour is ambiguous at saddle points. It is possible to resolve the ambiguity by using the average data value for the center of the cell to choose between different connections of the interpolated points. Here is another summary of the method showing options for the saddle : The central value is used to flip the index value before looking-up the cell geometry in the table, i.e. if the index is 0101=5 and the central value is below, then lookup index 10; similarly, if the index is 1010=10 and the central value is below, then lookup index 5.
Isoband
A similar algorithm can be created for filled contour bands within upper and lower threshold values. To build the index we compare the data values at the cell corners with the two contour threshold values. There are now 3 possibilities:
0 – corner data value below lower band level
1 – corner data value within band interval
2 – corner data value above upper level
The index will be ternary value built from these ternary digits, or trits. We build the index as before, by walking clockwise around the cell, appending each trit to the index, taking the most-significant-tritfrom the top left corner, and the least-significant-trit from the bottom left corner. There will now be 81 possibilities, rather than 16 for isolines. Each cell will be filled with 0, 1 or 2 polygonal fragments, each with 3–8 sides. The action for each cell is based on the category of the ternary index:
Empty – no fragments for index values 0 or 80.
Not Empty – walk around the cell adding corner positions that are withinthe band and interpolating along relevant edges; use the index to decide how to connect the list of points:
* Slope – build a single polygon fragment with 3–7 sides.
* Saddle – calculate the average value to help disambiguation; then use the index and the central value to build 1 or 2 polygonal fragments with a total of 6, 7 or 8 sides.
The case missing from the 6-sided saddles is for a central value that cannot occur. There is a valid case omitted from each 7-sided saddle, where the central value is dominated by a single extreme value. The resulting geometric structure would be too complex to fit the simple model of two convex fragments, so it is merged with the case where the central value is within the band. The linear interpolation in such cases will produce a plausible single heptagon.
Contouring triangle meshes
The same basic algorithm can be applied to triangular meshes, which consist of connected triangles with data assigned to the vertices. For example, a scattered set of data points could be connected with a Delaunay triangulation to allow the data field to be contoured. A triangular cell is always planar, because it is a 2-simplex. There is always a unique linear interpolant across a triangle and no possibility of an ambiguous saddle.
Isolines
The analysis for isolines over triangles is especially simple: there are 3 binary digits, so 8 possibilities:
Isobands
The analysis for isobands over triangles requires 3 ternary trits, so 27 possibilities:
Dimensions and spaces
The data space for the Marching Squares algorithm is 2D, because the vertices assigned a data value are connected to their neighbors in a 2D topological grid, but the spatial coordinates assigned to the vertices can be in 2D, 3D or higher dimensions. For example, a triangular mesh may represent a 2D data surface embedded in 3D space, where spatial positions of the vertices and interpolated points along a contour will all have 3 coordinates. Note that the case of squares is ambiguous again, because a quadrilateral embedded in 3-dimensional space is not necessarily planar, so there is a choice of geometrical interpolation scheme to draw the banded surfaces in 3D.
A naive implementation of Marching Squares that processes every cell independently will perform every linear interpolation twice or four times. Similarly, the output will contain 2 copies of the 2D vertices for disjoint lines or 4 copies for polygons. It is possible to reduce the computational overhead by caching the results of interpolation. For example, a single-threaded serial version would only need to cache interpolated results for one row of the input grid. It is also possible to reduce the size of the output by using indexed geometric primitives, i.e. create an array of 2D vertices and specify lines or polygons with short integer offsets into the array.
