Sort polygon's points for drawing

Question

I have a matrix (0 means nothing, 1 means terrain) that represents a level in my game. The matrix corresponds to a grid that my screen is broken up into, and indicates where my terrain goes.

My terrain is actually composed of 4 points in the corners of each block within the grid. When you have multiple blocks that are connected, I use a merge-cell algorithm that removes the duplicate points and any interior points. The result is that I end up with a list of points representing only the outer edges of the polygon.

To draw this polygon, I need the points to be in some sort of order (either clockwise or counter-clockwise) such that each point is followed by it's neighboring point. Obviously the first and last points need to be neighbors. Since this is all in a grid, I know the exact distance between neighboring points.

The problem is that I am having trouble coming up with an algorithm that allows me to "walk" around the edge of the polygon while putting the points in order. I believe there should be a way to utilize the fact that I have the matrix representing the geometry, meaning there is only 1 possible way to draw the polygon (even if it is concave).

I have tried several approaches using greedy-type algorithms, but can't seem to find a way to know, in every case, which direction I want to travel in. Given that any particular point can have up to 3 neighbors (the fourth isn't included because it is the "starting" point, meaning that I have already sorted it) I need a way of knowing which way to move.

Update

Another approach that I have been trying is to sort the points by their X (with tiebreaker of Y) which gives me the topmost/leftmost edge. It also guarantees that I am starting on an outer edge. However, I'm still struggling to find an algorithm that guarantees that I stay on the outside without crossing over.

Here is an example matrix:

0 0 0 0 0 0 0 0 0 0

0 0 0 1 1 1 0 0 0 0

0 0 0 0 0 1 0 0 0 0

0 0 0 0 0 1 1 1 0 0

Which corresponds to this (black dots represent my points):

Answer 1

First of all please consider that for a general matrix the output can be composed of more than one closed loop; for example boundaries of the matrix

form three distinct loops, one of them placed inside another.

To extract these loops the first step is to build a map of all "walls": you have a vertical wall each time the content of one cell is different from the next cell on the same row; you have instead an horizontal wall when the content is different from the same cell in the next row.

data = [[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ],
        [ 0, 1, 1, 1, 1, 0, 0, 0, 0, 0 ],
        [ 0, 1, 0, 0, 1, 0, 1, 1, 0, 0 ],
        [ 0, 1, 0, 0, 1, 0, 1, 1, 1, 0 ],
        [ 0, 1, 1, 1, 1, 0, 0, 1, 1, 0 ],
        [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ]]

rows = len(data)
cols = len(data[0])

walls = [[2*(data[r][c] != data[r][c+1]) + (data[r][c] != data[r+1][c])
          for c in range(cols-1)]
         for r in range(rows-1)]

In the example above I'm using two bits: 0x01 to mark horizontal walls and 0x02 to mark vertical walls. For a given (r, c) cell the walls are the right and bottom wall of the cell.

For simplicity I'm also assuming that the interesting areas are not touching the limits of the matrix; this can be solved by either adding extra rows and cols of zeros or by wrapping matrix access in a function that returns 0 for out-of-matrix virtual elements.

To build the list of boundaries you need to simply start from any point on a wall and move following walls, removing the walls from the map as you process them. When you cannot move any more a cycle has been completed (you're guaranteed to complete cycles because in a graph built in this way from a matrix of inside/outside flags the degree is guaranteed to be even in all vertices).

Filling all those cycles simultaneously using odd-even filling rules is also guaranteed to reproduce the original matrix.

In the code following I'm using r and c as row/col index and i and j instead to represent points on the boundary... for example for cell (r=3, c=2) the schema is:

where the red wall corresponds to bit 0x02 and the green wall to bit 0x01. The walls matrix has one row and one column less than the original data matrix because it's assumed that no walls can be present on last row or column.

result = []
for r in range(rows-1):
    for c in range(cols-1):
        if walls[r][c] & 1:
            i, j = r+1, c
            cycle = [(i, j)]
            while True:
                if i < rows-1 and walls[i][j-1] & 2:
                    ii, jj = i+1, j
                    walls[i][j-1] -= 2
                elif i > 0 and walls[i-1][j-1] & 2:
                    ii, jj = i-1, j
                    walls[i-1][j-1] -= 2
                elif j < cols-1 and walls[i-1][j] & 1:
                    ii, jj = i, j+1
                    walls[i-1][j] -= 1
                elif j > 0 and walls[i-1][j-1] & 1:
                    ii, jj = i, j-1
                    walls[i-1][j-1] -= 1
                else:
                    break
                i, j = ii, jj
                cycle.append((ii, jj))
            result.append(cycle)

Basically the code starts from a point on a boundary and the checks if it can move on a wall going up, down, left or right. When it cannot move any more a cycle has been completed and can be added to the final result.

The complexity of the algorithm is O(rows*cols), i.e. it's proportional to the input size and it's optimal (in big-O sense) because you cannot compute the result without at least reading the input. This is easy to see because the body of the while cannot be entered more times than the total number of walls in the map (at each iteration a wall is removed).

Edit

The algorithm can be modified to generate as output only simple cycles (i.e. paths in which each vertex is visited only once).

result = []
index = [[-1] * cols for x in range(rows)]
for r in range(rows-1):
    for c in range(cols-1):
        if walls[r][c] & 1:
            i, j = r+1, c
            cycle = [(i, j)]
            index[i][j] = 0
            while True:
                if i > 0 and walls[i-1][j-1] & 2:
                    ii, jj = i-1, j
                    walls[i-1][j-1] -= 2
                elif j > 0 and walls[i-1][j-1] & 1:
                    ii, jj = i, j-1
                    walls[i-1][j-1] -= 1
                elif i < rows-1 and walls[i][j-1] & 2:
                    ii, jj = i+1, j
                    walls[i][j-1] -= 2
                elif j < cols-1 and walls[i-1][j] & 1:
                    ii, jj = i, j+1
                    walls[i-1][j] -= 1
                else:
                    break
                i, j = ii, jj
                cycle.append((ii, jj))
                ix = index[i][j]
                if ix >= 0:
                    # closed a loop
                    result.append(cycle[ix:])
                    for i_, j_ in cycle[ix:]:
                        index[i_][j_] = -1
                    cycle = cycle[:ix+1]
                index[i][j] = len(cycle)-1

This is implemented by adding to the output a separate cycle once the same vertex is met twice in the processing (the index table stores for a given i,j point the 0-based index in the current cycle being built).