Determining orthogonal polygons in a 3D integer lattice The last couple of days I've been trying to get a simple algorithm for uniquely determining orthogonal polygons in 3D integer space, having already determined coplanar points.
By the time we got here we know each point is unique and there are no collinear points.
For example, given the coplanar points pts (lists or tuples)
    pts = [
        [0, 0, 0], [4, 0, 0],
        [1, 1, 0], [2, 1, 0], [3, 1, 0], [4, 1, 0],
        [1, 2, 0], [2, 2, 0], [3, 2, 0], [4, 2, 0],
        [0, 3, 0], [4, 3, 0]
    ]

# Generate orthogonal polygons
# 3 +-----------+
#   |           |
# 2 |  +__+  +--+
#   |  |  |  |
# 1 |  +--+  +__+
#   |           |
# 0 +-----------+
#   0  1  2  3  4

#    result = [
#       [
#         [0, 0, 0], [0, 3, 0], [4, 3, 0], [4, 2, 0], 
#         [3, 2, 0], [3, 1, 0], [4, 1, 0], [4, 0, 0]
#       ], [
#         [1, 1, 0], [1, 2, 0], [2, 2, 0], [2, 1, 0]
#       ]
#    ]


I found the basis of an algorithm for this in a paper as authored by Joseph O'Rourke
The following seems to work for the few tests I set up (however it doesn't deal with orientation).

def polygonise(pts: list) -&gt; list:
    # axes gives us the varying columns (we can ignore z for this)
    axes = {i:set(c) for i,c in enumerate(zip(*pts)) if len(set(c))!=1}
    a, b = axes.keys()

    # construct dicts of axis-columns 
    pts_by_axis = {axis: {col: [] for col in cols} for axis,cols in axes.items()}
    for pt in pts:
        for axis in axes:
            pts_by_axis[axis][pt[axis]].append(pt)

    # each axis-column (sorted by the other dimension)
    # will now be paired off into edges, and constructed into axis-edge-dicts
    edges = {}
    for axis in axes:
        other = a if axis != a else b
        for col,pt_list in pts_by_axis[axis].items():
            pt_list.sort(key=lambda p:p[other])
            for i in range(0,len(pt_list),2):
                v1 = tuple(pt_list[i])
                v2 = tuple(pt_list[i + 1])
                edges[tuple((axis,v1))] = other, v2
                edges[tuple((axis,v2))] = other, v1


    # now derive the polygons, starting at any point.
    axis, edge = next(iter(edges.keys()))
    result, poly = [], [edge]

    while edges:
        if (axis, edge) not in edges:
            # This polygon is now closed. keep it, and move to the next one.
            result.append(poly[:-1])  # if wanted looped, don't slice
            poly = []
            axis, edge = next(iter(edges.keys()))
        else:
            # Get the next edge. change the axis (edge is on a different axis)
            axis, new_edge = edges.pop((axis,edge))
            # delete any old reversed edge.
            del edges[(axis, edge)]
            edge = new_edge

        # add the edge to the existing polygon
        poly.append(edge)
    # capture the final polygon.
    result.append(poly[:-1]) # if wanted looped, don't slice.
    return result



if __name__ == &quot;__main__&quot;:
    from mpl_toolkits.mplot3d.art3d import Poly3DCollection
    import matplotlib.pyplot as plt
    from random import shuffle
    pts = [
        (0, 0, 0), (4, 0, 0),
        (1, 1, 0), (2, 1, 0), (3, 1, 0), (4, 1, 0),
        (1, 2, 0), (2, 2, 0), (3, 2, 0), (4, 2, 0),
        (0, 3, 0), (4, 3, 0)
    ]
    # shuffle just to make it different.
    shuffle(pts)
    result = polygonise(pts)
    # all the below is just to show the result
    fig = plt.figure()
    subplot = fig.add_subplot(111, projection='3d')
    subplot.axis('off')
    subplot.scatter([0.6, 3.1], [0.7, 3.3], [0, 0], alpha=0.0)
    subplot.add_collection3d(Poly3DCollection(result, edgecolors='k', facecolors='b', linewidths=1, alpha=0.2))
    fig.tight_layout()
    plt.show()



I think that the loops could be easier, and I'm not convinced that I need to store edges twice, (which means I have to delete the unused one).
Any ideas for tidying the polygonise() method would be welcome.
J. O'Rourke, &quot;Uniqueness of orthogonal connect-the-dots&quot;, Computational Morphology, G.T. Toussaint (editor), Elsevier Science Publishers, B.V.(North-Holland), 1988, 99-104. found here
Type Hints
+1 for using type hints.  But they could be a little more informative.
def polygonise(pts: list) -&gt; list:
Ok.  The function takes a list of pts, but what are pts?  It returns a list, but a list of what?
from typing import Sequence

Point = Sequence[int]
Polygon = list[Point]

def polygonise(points: list[Point]) -&gt; list[Polygon]:
    ...

Now we're giving the reader a lot of information.  A point is a sequence of integers.  Could be a tuple, could be a list.  Not important.  Instead of pts, we're a bit more explicit and say points.  Finally, the function returns a specific kind of list: a Polygon list.
You can help the reader out even more by including a nice &quot;&quot;&quot;docstring&quot;&quot;&quot; for the function.
Axes
axes = {i:set(c) for i,c in enumerate(zip(*pts)) if len(set(c))!=1}
I really, really hate this statement.  It takes the list of points, and explodes each point as an argument for zip.  Then, it takes the first coordinate of each of those points, and returns it as a tuple.  enumerate adds an index to this, to keep track of which coordinate axis we are processing.  Then, we assign those to things to i, c ... really opaque variable names.  axis, coordinates might be clearer.  Then we turn the coordinates into a set, determine the len of the set, and if that is not 1, we turn the coordinates into a set again, for construction the dictionary.
Yikes!  That's a lot to unpack.
Then we take the axes, and assign it into to variables a and b, hoping there are exactly 2 axis.
Perhaps we should add some error checking.
    num_coord = len(points[0])
    if any(len(pt) != num_coord for pt in points):
        raise ValueError(&quot;All points should have the same number of coordinates&quot;)

    axes = ...

    if len(axes) != 2:
        raise ValueError(&quot;All coordinates must be in an orthogonal 2D plane&quot;)

    a, b = axes.keys()

That gives me a lot more confidence in the code.
I still don't like the axes = ... statement.  I would re-write it, and maybe use the walrus operator (:=) to avoid creating the set twice, but I'm about to significantly change it.
Useless code
edges[tuple((axis,v1))] = other, v2
Here, (axis, v1) is creating a tuple, to pass to the tuple(...) function.  That seems ... odd.  You could simplify it to:
edges[axis, v1] = other, v2
Alternate implementation
Axes
Instead of a dictionary, let's just extract a list of axis which have varying coordinate data.  A simple comparison of each coordinate of all points with the coordinate of the first point, for each coordinate, will tell us if that axis is important or not.
    axes = [axis for axis, first in enumerate(points[0])
            if any(pt[axis] != first for pt in points)]

    if len(axes) != 2:
        raise ValueError(&quot;All coordinates must be in an orthogonal 2D plane&quot;)

    a, b = axes

This is a much simpler concept.  No sets are being created (twice).  No dictionaries.
Sorted Edges
The pts_by_axis is very complex, a dictionary of dictionaries of arrays.
Let's rework this, by simply sorting the points list.  We'll sort it two ways, first by the (a, b) coordinate, then by the (b, a) coordinate.  We'll use a helper function to generate the sort keys:
def point_key(first: int, second: int):
    def key(pt):
        return pt[first], pt[second]
    return key

Now we can say sorted(points, key=point_key(a, b)), and we'll have the points sorted in one direction.  sorted(points, key=point_key(b, a)) will sort it in the other direction.
Since the points always come in pairs, to form edges, let's add another helper function, to extract pairs of points from the sorted array.
def in_pairs(iterable):
    it = iter(iterable)
    return zip(it, it)

With these, we can build our sorted edge structure.  But instead of keying the edges dictionary with a tuple to encoding the axis of the edge, let's just use two dictionaries:
    a_edge = {}
    for v1, v2 in in_pairs(sorted(points, key=point_key(a, b))):
        a_edge[v1] = v2
        a_edge[v2] = v1

    b_edge = {}
    for v1, v2 in in_pairs(sorted(points, key=point_key(b, a))):
        b_edge[v1] = v2
        b_edge[v2] = v1

The loop at the end does some check to close a polygon and start the next one.  Really, you have two loops.  Loop while there are still edges.  Inside that loop, loop until you've made a polygon.  With that structure, the code reads much better:
    result = []
    while a_edge:
        v1 = next(iter(a_edge.keys()))
        poly = []

        while v1 in a_edge:
            v2 = a_edge.pop(v1)
            v3 = b_edge.pop(v2)
            a_edge.pop(v2)
            b_edge.pop(v1)
            poly.extend((v1, v2))
            v1 = v3

        result.append(poly)
        
    return result

Note that I use poly = [] instead of poly = [v1], which means no slice is required at the end to remove the extra point.
Also, I'm popping v2 out of the a_edge dictionary, not a edge, since what is popped out is just a vertex.  I find that much clearer.