连接点网络中的多边形

Question

给定一个 2D 点数组 (#pts x 2) 和一个由哪些点连接到哪个点的数组（#bonds x 2 int 数组，索引为 pts），我如何有效地返回一个由债券？

可能存在不闭合多边形的“悬空”键（如下图左上角所示），这些应被忽略。

这是一个例子：

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

对于未来的读者，弗兰克的建议在 numpy 中的大部分实现如下。边界的提取基本上遵循与绕多边形行走相同的算法，不同之处在于使用最小角度键，而不是最大值。

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

Answer 2

我不能告诉你如何用 numpy 来实现它，但这里有一个可能的算法的概要：

1. 为每个点添加一个附加键列表。
2. 删除仅连接一个键的点，同时删除此键（这些是悬空键）
3. 为每个键附加两个布尔标记，指示键是否已在两个可能方向之一上添加到多边形。每个键只能用于两个多边形。最初将所有标记设置为 false。
4. 选择任意初始点并重复以下步骤，直到在两个方向上都使用了所有键：
5. 选择未使用的键（在相应方向上）。这是多边形的第一条边。在连接到所选端点的键中，选择角度最小的键，例如逆时针方向。将此添加到多边形并继续，直到返回初始点。

此算法还将生成一个包含网络所有外部键的大多边形。我想你会找到一种方法来识别并删除它。