How to order points for Matplotlib's fill function?

Question

A polygon is represented by two lists. List xs contains x coordinates in decreasing order. For each x value, the corresponding element in list regions is a list of slices. Each slice indicates a y range that belongs to the interior of a polygon. For each x value, the slices are sorted in increasing order.

I want to draw the filled polygon using Matplotlib's fill function. That function needs the points to be ordered so that going from point to point describes the contour of the polygon. How can I order the points in my data set to obtain the right polygon?

Here is an example of the kind of data I have.

xs = range(10, 0, -1)
regions = [[slice(0, 3)],
           [slice(0, 4), slice(5.2, 5.8)],
           [slice(1, 5), slice(5.4, 5.8)],
           [slice(1.3, 6)],
           [slice(1.8, 6)],
           [slice(1.9, 3), slice(3.1, 6.1)],
           [slice(2, 2.9), slice(3.2, 5), slice(6, 6.2)],
           [slice(2.1, 2.7), slice(3.4, 5.1), slice(5.9, 6.3)],
           [slice(3.5, 5.2), slice(5.7, 6.4)],
           [slice(3.8, 5.3), slice(5.8, 6.1)]]

The correct way to order the points is as follows

xx = [10, 9, 8, 7, 6, 5, 4, 3, 3, 4, 5, 5, 4, 3, 2, 1,
   1, 2, 3, 4, 4, 3, 2, 1, 1, 2, 3, 4, 5, 6, 7, 8,
   9, 9, 8, 8, 9, 10]
yy = [0, 0, 1, 1.3, 1.8, 1.9, 2, 2.1, 2.7, 2.9, 3, 3.1,
      3.2, 3.4, 3.5, 3.8, 5.3, 5.2, 5.1, 5, 6, 5.9, 5.7,
      5.8, 6.1, 6.4, 6.3, 6.2, 6.3, 6, 6, 5.8, 5.8, 5.2,
      5.4, 5, 4, 3]

and the figure should look as such

So far, I have tried defining turning points, i.e., x values where the contour changes direction. This can be done with numpy.diff applied to an array containing the number of slices for each x value. If the difference is not zero, that x value is a turning point. This can probably be used to figure out the next point. The difficulty is to figure out which slice is next and whether to use the start or the end of the slice.

This question is similar, but in my case the polygons have a much more complicated shape. On the other hand, I do have information about what is inside the polygon.

Edit

I just realized that the problem as stated does not always have a unique solution. For instance, the set of points in the example above also admits the solution

xx = [10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 1, 2, 3, 4, 5, 5,
       4, 3, 3, 4, 4, 3, 2, 1, 1, 2, 3, 4, 5, 6, 7, 8,
       9, 9, 8, 8, 9, 10]
yy = [0, 0, 1, 1.3, 1.8, 1.9, 2, 2.1, 3.5, 3.8, 5.3, 5.2,
      2.7, 2.9, 3, 3.1, 3.2, 3.4, 5.1, 5, 6, 5.9, 5.7,
      5.8, 6.1, 6.4, 6.3, 6.2, 6.1, 6, 6, 5.8, 5.8, 5.2,
      5.4, 5, 4, 3]

The solution I am looking for is the one that minimizes the absolute value of each edge's slope (apart from the vertical segments). In the example above, this corresponds to the first solution.

Answer 1

Here is function that does the job. It is a generator that yields all polygons found (in case there is more than one connected polygon).

First, the regions have to be rewritten as tuples:

regions = [[(sl.start, sl.stop) for sl in sllist] for sllist in regions]

Ordering is done by calling the generator (in a loop or otherwise). The function yields a polygon as a list of (x, y) points. Obtaining the separate xx and yy lists is easy:

for polygon in order_points(xs, regions):
    xx, yy = zip(*polygon)

The function itself is commented in detail.

def order_points(x, regions):
    """Given two lists of the same length, one of x values and one of regions,
    determine how to order the points so that they describe a polygon that
    contains the interior of all regions.

    This function is a generator that yields disjoint polygons.

    Parameters
    ==========

    x: list
       x values at which the regions are defined

    regions: list of lists of tuples
       Each list contains information for one x value.  The sublist contains
       so-called regions that are described by tuples corresponding to the
       bottom and the top of the region.  The y values between the bottom and
       the top of each region are in the interior of a polygon.

       Each list of tuples should be sorted by increasing y value.

    Yields
    ======

    polygon: list of tuples
       Ordered list of (x, y) coordinated of the points on the polygon.

    """
    # Make sure the x values are in ascending order.
    xvals, yregions = zip(*sorted(zip(x, regions)))

    # Direction is -1 when going toward low x, 1 when going toward high x.
    direction = 1
    # Indicate whether the inside of the polygon is below, 0, or above, 1, the
    # current point.
    inside = 1

    # List all possible combinations of x index, region index.
    tovisit = [(pos, rpos) for pos in range(len(xvals))
                            for rpos in range(len(yregions[pos]))]
    pos, rpos = tovisit.pop(0)
    ycur = yregions[pos][rpos][0]

    # Keep track of already visited points.
    visited = set()
    # Keep track of current polygon.
    polygon = []
    while True:
        # Keep track of ordered points.
        xcur = xvals[pos]
        polygon.append((xcur, ycur))
        visited.add((xcur, ycur))

        # Find the minimum vertical distance between the current position and
        # the following points: points at the next (as specified by direction)
        # x value, and points at the current x value

        # For points at the next x, if the polygon is currently above, inside =
        # 1, the points considered are at the bottom of a region, i.e., at
        # index 0 of the tuple.
        next_pos = pos + direction
        if next_pos < 0 or next_pos >= len(xvals):
            next_pos = pos
        distance = -1
        for ri, region in enumerate(yregions[pos]):
            if (xcur, region[inside]) not in visited:
                d = abs(ycur - region[inside])
                if d < distance or distance == -1:
                    distance = d
                    move = ('vertical', (pos, ri, inside))

        # For points at the same x, if the polygon is currently above, the
        # points considered are at the top of a region, i.e., at index 1 of the
        # tuple.
        for ri, region in enumerate(yregions[next_pos]):
            polypos = (not inside)
            if (xvals[next_pos], region[polypos]) not in visited:
                d = abs(ycur - region[polypos])
                if d < distance or distance == -1:
                    distance = d
                    move = ('next', (next_pos, ri, polypos))

        # If no suitable next point is found, the polygon is complete. Yield
        # the polygon and try to find a separate polygon.
        if distance < 0:
            yield polygon
            polygon = []
            direction = 10  # dummy value to detect if no next point is found
            while tovisit:
                pos, slpos = tovisit.pop(0)
                ycur = yregions[pos][slpos][0]
                if (xvals[pos], ycur) not in visited:
                    direction = 1
                    inside = 1
                    break
            if direction == 10:
                break
            continue

        # Make the move.
        if move[0] == 'vertical':
            # Vertical move changes the direction (unless it is the very first
            # move) and toggles inside/outside.
            if len(polygon) > 1:
                direction *= -1
            inside = int(not inside)
        pos, rpos, end = move[1]
        ycur = yregions[pos][rpos][end]