Get the coordinates from ContourPlot and RegionPlot

Question

How do I get the coordinates from a contour plot I've done in Mathematica? For example, I have a two-variable function f[x, y], for which I can make a contour plot:

    contour = ContourPlot[f[x, y] == 1, {x, -1, 1}, {y, -1, 1}];

I can access a nested list, containing lists of coordinates, from an exported file using

    Export["output.m", contour, "TEXT"] 

but the nested list is inside the function Graphics.

I would like to export only the corresponding nested list. What is a straightforward way to do that? Is there a way of doing that by manipulating directly the variable contour?

Edit

@Jens solution was

    Cases[Normal@contour, Line[x_] :> x, Infinity]

Is there a way to generalize it? The issue is that Line won't give us the points inside a region. E.g., I would like to get the coordinates corresponding to a region plot:

    region = RegionPlot[f[x, y] > 1, {x, -1, 1}, {y, -1, 1}];

Answer 1

Since the plot usually is a GraphicsComplex, the extraction is easiest if you first convert using Normal:

contour = ContourPlot[x^2 + y^2 == 1, {x, -1, 1}, {y, -1, 1}];

Cases[Normal@contour, Line[x_] :> x, Infinity]

This produces a list that shows the coordinates in the order they were drawn.

Explanation:

The contours in the plot are drawn using the Line command, which takes lists of points (or lists of lists of points). Usually, these points are all collected at the front of the ContourPlot output in the form of a GraphicsComplex, such that each point can later be addressed by using an index from within the Line commands. By applying Normal, these indexed points are moved to where they are actually used in the drawing part of the output. Normal@contour is the same as Normal[contour].

After that is done, we can look for all Line commands and find the coordinates in sequential order inside of them. This is done by using Cases, which selects parts of the expression that match a pattern. The pattern is specified here as Line[x_] where x_ is a "dummy variable" that gets defined whenever a Line was found, by replacing it with the contents of the line. The final step is to tell Cases that when it does find a value for x, to just output that without the wrapper Line.

This search is done throughout the whole plot, which is indicated by the Infinity level specification.

Update for RegionPlot

Extracting points from a plot using Cases can be generalized to situations where the points aren't inside a Line. In RegionPlot, for example, you may want to extract the points that form the mesh with which the region is filled. This filling is done by a polygonal tesselation, so we have to simply replace Line with Polygon:

region = Normal@RegionPlot[x^2 + y^2 <= 1, {x, -1, 1}, {y, -1, 1}];
pts = DeleteDuplicates@
   Flatten[Cases[region, Polygon[x_] :> x, Infinity], 1];

Graphics[Point[pts]]

Here I had to add two more steps before being able to make the plot: first, Flatten is used to remove all nested levels except the ones grouping the coordinate tuples for each point. Then, I added DeleteDuplicates just to remove any shared vertices between polygons, so one and the same point isn't re-drawn redundantly.