How does Matlab calculate contour lines?

Question

What is the algorithm that Matlab uses to generate contour lines? In other words, how does it transform level data on a grid into a set of lines?

What I would like is:

• the local criterion to obtain points that lie on the contour?
• the global procedure to capture all contour lines?

I don't need detailed specifics about the underlying code, but the general principles would be helpful for me to interpret the output. I use contour (and derivatives) in my research, and want to get a sense of the numerical errors that are introduced in this step.

It looks like a very simple question, but I couldn't find an explanation in Matlab's documentation, nor found anything on SO or elsewhere on the web. My apologies if it turns about to be easy to find after all.

Answer 1

From MATLAB® - Graphics - R2012a, from page 5-73 to page 5-76:

The Contouring Algorithm

The contourc function calculates the contour matrix for the other contour functions. It is a low-level function that is not called from the command line. The contouring algorithm first determines which contour levels to draw. If you specified the input vector v, the elements of v are the contour level values, and length(v) determines the number of contour levels generated. If you do not specify v, the algorithm chooses no more than 20 contour levels that are divisible by 2 or 5.

The height matrix Z has associated X and Y matrices that locate each value in Z at the intersection of a row and a column, or these matrices are inferred when they are unspecified. The row and column widths can vary, but typically they are constant (i.e., Z is a regular grid). Before calling contourc to interpolate contours, contourf pads the height matrix with an extra row or column on each edge. It assigns z-values to the added grid cells that are well below the minimum value of the matrix. The padded values enable contours to close at the matrix boundary so that they can be filled with color. When contourc creates the contour matrix, it replaces the x,y coordinates containing the low z-values with NaNs to prevent contour lines that pass along matrix edges from being displayed. This is why contour matrices returned by contourf sometimes contain NaN values. Set the current level, c, equal to the lowest contour level to be plotted within the range [min(Z) max(Z)]. The contouring algorithm checks each edge of every square in the grid to see if c is between the two z values for the edge points. If so, a contour at that level crosses the edge, and a linear interpolation is performed:

t=(c-Z0)/(Z1-Z0)

Z0 is the z value at one edge point, and Z1 is the z value at the other edge point.

Start indexing a new contour line (i=1) for level c by interpolating x and y:

cx(i) = X0+t*(X1-X0)
cy(i) = Y0+t*(Y1-Y0)

Walk around the edges of the square just entered; the contour exits at the next edge with z values that bracket c. Increment i, compute t for the edge, and then compute cx(i) and cy(i), as above. Mark the square as having been visited. Keep checking the edges of each square entered to determine the exit edge until the line(cx,cy) closes on its initial point or exits the grid. If the square being entered is already marked, the contour line closes there. Copy cx, cy, c, and i to the contour line data structure (the matrix returned by contouring functions, described shortly).

Reinitialize cx, cy, and i. Move to an unmarked square and test its edges for intersections; when you find one at level c, repeat the preceding operations. Any number of contour lines can exist for a given level. Clear all the markers, increment the contour level, and repeat until c exceeds max(Z). Extra logic is needed for squares where a contour passes through all four edges (saddle points) to determine which pairs of edges to connect. contour, contour3, and contourf return a two-row matrix that specifies all the contour lines:

C = [ value1 xdata(1) xdata(2)...
        numv ydata(1) ydata(2)...]

The first row of the column that begins each definition of a contour line contains the contour value, as specified by v and used by clabel. Beneath that value is the number of (x,y) vertices in the contour line. Remaining columns contain the data for the (x,y) pairs. For example, the contour matrix calculated by C = contour(peaks(3)) is as follows.

The circled values begin each definition of a contour line.

Answer 2

You can read about Marching Squares. (https://en.wikipedia.org/wiki/Marching_squares)

Usually, linear interpolation is performed along the grid cell edges, giving contour points that you link to form polylines.

When the cells are too coarse, a preliminary interpolation, for instance bicubic, can be made to refine the mesh.