Determine order of 2-D convex hull points based on the neighborhood structure

cvx2d.py

import numpy as np
from scipy.spatial import ConvexHull
from scipy.sparse import coo_matrix, csgraph

def ordered_hull_idx_2d(hull):
    n = hull.simplices.shape[0]
    
    # determine order of edges in the convex hull
    v = coo_matrix((np.ones(2*n), (np.repeat(np.arange(n), 2), hull.neighbors.ravel())))
    facet_order = csgraph.depth_first_order(v, 0, return_predecessors=False)
    facet_vidx = hull.simplices[facet_order]
    
    # pick one vertex for each edge, based on which direction the walk went
    m = hull.neighbors[facet_order][:-1] == facet_order[1:,None]
    i = np.arange(n)
    j = np.r_[np.where(m)[1], 0] 
    
    ordered_vertex_idx = facet_vidx[i, j]
    
    # sanity check
    assert np.all(np.unique(ordered_vertex_idx) == np.unique(hull.simplices.ravel()))

    return ordered_vertex_idx

def main():
    points = np.random.rand(1000, 2)
    hull = ConvexHull(points)
    hull_idx = ordered_hull_idx_2d(hull)

    # Close the polygon
    hull_idx = np.r_[hull_idx, hull_idx[0]]

    # Plot the points
    hull_pts = points[hull_idx]

    import matplotlib.pyplot as plt
    plt.clf()
    plt.plot(points[:,0], points[:,1], '.',
             hull_pts[:,0], hull_pts[:,1], '-')
    plt.xlim(-0.5, 1.5)
    plt.ylim(-0.5, 1.5)
    plt.show()

    return hull

if __name__ == "__main__":
    main()

cvx2dcy.pyx

# The cython implementation is simpler, we can just walk around

import numpy as np

def ordered_hull_idx_2d(hull):
    cdef int[:,:] simplices = hull.simplices
    cdef int[:,:] neighbors = hull.neighbors
    cdef int[:] ordered
    cdef int j, n, ifacet, ifacet2, ipos

    assert simplices.shape[1] == 2

    n = simplices.shape[0]
    ordered = np.zeros(n, dtype=np.intc)

    ifacet = 0
    ipos = 0
    for j in range(n):
        # record current vertex idx
        ordered[j] = simplices[ifacet, ipos]

        # jump to the neighboring facet sharing the just
        # recorded vertex
        ifacet2 = neighbors[ifacet, 1-ipos]

        # the indices in the neighboring facet may be flipped,
        # so check which side is which, and pick the different
        # vertex
        if simplices[ifacet2, 0] == simplices[ifacet, ipos]:
            ipos = 1
        else:
            ipos = 0
        ifacet = ifacet2

    return np.asarray(ordered)