Polygons from network of connected points

Question

Given an array of 2D points (#pts x 2) and an array of which points are connected to which (#bonds x 2 int array with indices of pts), how can I efficiently return an array of polygons formed from the bonds?

There can be 'dangling' bonds (like in the top left of the image below) that don't close a polygon, and these should be ignored.

Here's an example:

import numpy as np
xy = np.array([[2.72,-2.976], [2.182,-3.40207],
[-3.923,-3.463], [2.1130,4.5460], [2.3024,3.4900], [.96979,-.368],
[-2.632,3.7555], [-.5086,.06170], [.23409,-.6588], [.20225,-.9540],
[-.5267,-1.981], [-2.190,1.4710], [-4.341,3.2331], [-3.318,3.2654],
[.58510,4.1406], [.74331,2.9556], [.39622,3.6160], [-.8943,1.0643],
[-1.624,1.5259], [-1.414,3.5908], [-1.321,3.6770], [1.6148,1.0070],
[.76172,2.4627], [.76935,2.4838], [3.0322,-2.124], [1.9273,-.5527],
[-2.350,-.8412], [-3.053,-2.697], [-1.945,-2.795], [-1.905,-2.767],
[-1.904,-2.765], [-3.546,1.3208], [-2.513,1.3117], [-2.953,-.5855],
[-4.368,-.9650]])

BL= np.array([[22,23], [28,29], [8,9],
[12,31], [18,19], [31,32], [3,14],
[32,33], [24,25], [10,30], [15,23],
[5,25],  [12,13], [0,24],  [27,28],
[15,16], [5,8],   [0,1],   [11,18],
[2,27],  [11,13], [33,34], [26,33],
[29,30], [7,17],  [9,10],  [26,30],
[17,22], [5,21],  [19,20], [17,18],
[14,16], [7,26],  [21,22], [3,4],
[4,15],  [11,32], [6,19],  [6,13],
[16,20], [27,34], [7,8],   [1,9]])

Answer 1

I can't tell you how to implement it with numpy, but here's an outline of a possible algorithm:

1. Add a list of attached bonds to each point.
2. Remove the points that have only one bond attached, remove this bond as well (these are the dangling bonds)
3. Attach two boolean markers to each bond, indicating if the bond has already been added to a polygon in one of the two possible directions. Each bond can only be used in two polygons. Initially set all markers to false.
4. Select any initial point and repeat the following step until all bonds have been used in both directions:
5. Select a bond that has not been used (in the respective direction). This is the first edge of the polygon. Of the bonds attached to the end point of the selected one, choose the one with minimal angle in e.g. counter-clockwise direction. Add this to the polygon and continue until you return to the initial point.

This algorithm will also produce a large polygon containing all the outer bonds of the network. I guess you will find a way to recognize this one and remove it.

Answer 2

For future readers, the bulk of the implementation of Frank's suggestion in numpy is below. The extraction of the boundary follows essentially the same algorithm as walking around a polygon, except using the minimum angle bond, rather than the max.

def extract_polygons_lattice(xy, BL, NL, KL):
    ''' Extract polygons from a lattice of points.

    Parameters
    ----------
    xy : NP x 2 float array
        points living on vertices of dual to triangulation
    BL : Nbonds x 2 int array
        Each row is a bond and contains indices of connected points
    NL : NP x NN int array
        Neighbor list. The ith row has neighbors of the ith particle, padded with zeros
    KL : NP x NN int array
        Connectivity list. The ith row has ones where ith particle is connected to NL[i,j]

    Returns
    ----------
    polygons : list
        list of lists of indices of each polygon
    PPC : list
        list of patches for patch collection
    '''
    NP = len(xy)
    NN = np.shape(KL)[1]
    # Remove dangling bonds
    # dangling bonds have one particle with only one neighbor
    finished_dangles = False
    while not finished_dangles:
        dangles = np.where([ np.count_nonzero(row)==1 for row in KL])[0]
        if len(dangles) >0:
            # Make sorted bond list of dangling bonds
            dpair = np.sort(np.array([ [d0, NL[d0,np.where(KL[d0]!=0)[0]] ]  for d0 in dangles ]), axis=1)
            # Remove those bonds from BL
            BL = setdiff2d(BL,dpair.astype(BL.dtype))
            print 'dpair = ', dpair
            print 'ending BL = ', BL
            NL, KL = BL2NLandKL(BL,NP=NP,NN=NN)
        else:
            finished_dangles = True


    # bond markers for counterclockwise, clockwise
    used = np.zeros((len(BL),2), dtype = bool)
    polygons = []
    finished = False

    while (not finished) and len(polygons)<20:
        # Check if all bond markers are used in order A-->B
        todoAB = np.where(~used[:,0])[0]
        if len(todoAB) > 0:
            bond = BL[todoAB[0]]

            # bb will be list of polygon indices
            # Start with orientation going from bond[0] to bond[1]
            nxt = bond[1]
            bb = [ bond[0], nxt ]
            dmyi = 1

            # as long as we haven't completed the full outer polygon, add next index
            while nxt != bond[0]:
                n_tmp = NL[ nxt, np.argwhere(KL[nxt]).ravel()]
                # Exclude previous boundary particle from the neighbors array, unless its the only one
                # (It cannot be the only one, if we removed dangling bonds)
                if len(n_tmp) == 1:
                    '''The bond is a lone bond, not part of a triangle.'''
                    neighbors = n_tmp
                else:
                    neighbors = np.delete(n_tmp, np.where(n_tmp == bb[dmyi-1])[0])

                angles = np.mod( np.arctan2(xy[neighbors,1]-xy[nxt,1],xy[neighbors,0]-xy[nxt,0]).ravel() \
                        - np.arctan2( xy[bb[dmyi-1],1]-xy[nxt,1], xy[bb[dmyi-1],0]-xy[nxt,0] ).ravel(), 2*np.pi)
                nxt = neighbors[angles == max(angles)][0]
                bb.append( nxt )


                # Now mark the current bond as used
                thisbond = [bb[dmyi-1], bb[dmyi]]
                # Get index of used matching thisbond
                mark_used = np.where((BL == thisbond).all(axis=1))
                if len(mark_used)>0:
                    #print 'marking bond [', thisbond, '] as used'
                    used[mark_used,0] = True
                else:
                    # Used this bond in reverse order
                    used[mark_used,1] = True

                dmyi += 1

            polygons.append(bb) 

        else:
            # Check for remaining bonds unused in reverse order (B-->A)
            todoBA = np.where(~used[:,1])[0]
            if len(todoBA) >0:
                bond = BL[todoBA[0]]

                # bb will be list of polygon indices
                # Start with orientation going from bond[0] to bond[1]
                nxt = bond[0]
                bb = [ bond[1], nxt ]
                dmyi = 1

                # as long as we haven't completed the full outer polygon, add nextIND
                while nxt != bond[1]:
                    n_tmp = NL[ nxt, np.argwhere(KL[nxt]).ravel()]
                    # Exclude previous boundary particle from the neighbors array, unless its the only one
                    # (It cannot be the only one, if we removed dangling bonds)
                    if len(n_tmp) == 1:
                        '''The bond is a lone bond, not part of a triangle.'''
                        neighbors = n_tmp
                    else:
                        neighbors = np.delete(n_tmp, np.where(n_tmp == bb[dmyi-1])[0])

                    angles = np.mod( np.arctan2(xy[neighbors,1]-xy[nxt,1],xy[neighbors,0]-xy[nxt,0]).ravel() \
                            - np.arctan2( xy[bb[dmyi-1],1]-xy[nxt,1], xy[bb[dmyi-1],0]-xy[nxt,0] ).ravel(), 2*np.pi)
                    nxt = neighbors[angles == max(angles)][0]
                    bb.append( nxt )


                    # Now mark the current bond as used --> note the inversion of the bond order to match BL
                    thisbond = [bb[dmyi], bb[dmyi-1]]
                    # Get index of used matching [bb[dmyi-1],nxt]
                    mark_used = np.where((BL == thisbond).all(axis=1))
                    if len(mark_used)>0:
                        used[mark_used,1] = True

                    dmyi += 1

                polygons.append(bb) 
            else:
                # All bonds have been accounted for
                finished = True


    # Check for duplicates (up to cyclic permutations) in polygons
    # Note that we need to ignore the last element of each polygon (which is also starting pt) 
    keep = np.ones(len(polygons),dtype=bool)
    for ii in range(len(polygons)):
        polyg = polygons[ii]
        for p2 in polygons[ii+1:]:
            if is_cyclic_permutation(polyg[:-1],p2[:-1]):
                keep[ii] = False

    polygons = [polygons[i] for i in np.where(keep)[0]]

    # Remove the polygon which is the entire lattice boundary, except dangling bonds
    boundary = extract_boundary_from_NL(xy,NL,KL)
    print 'boundary = ', boundary
    keep = np.ones(len(polygons),dtype=bool)
    for ii in range(len(polygons)):
        polyg = polygons[ii]
        if is_cyclic_permutation(polyg[:-1],boundary.tolist()):
            keep[ii] = False
        elif is_cyclic_permutation(polyg[:-1],boundary[::-1].tolist()):
            keep[ii] = False

    polygons = [polygons[i] for i in np.where(keep)[0]]


    # Prepare a polygon patch collection
    PPC = []  
    for polyINDs in polygons:
        pp = Path(xy[polyINDs],closed=True)
        ppp = patches.PathPatch(pp, lw=2)
        PPC.append(ppp)


    return polygons, PPC