來自連接點網絡的多邊形

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的實現它，但這裏有一個可能的算法的輪廓：

1. 添加附名單債券到每個點。
2. 刪除只有一個債券連接的點，也刪除這個債券（這些是懸掛債券）
3. 附加兩個布爾標記到每個債券，指示債券是否已添加到多邊形中的一個兩個可能的方向。每個債券只能用於兩個多邊形。最初將所有標記設置爲false。
4. 選擇任意初始點，並重復下列步驟，直到所有的鍵都在兩個方向上被使用：
5. 選擇未被使用（在各個方向）的鍵。這是多邊形的第一條邊。在連接到選定終點的債券中，選擇例如具有最小角度的債券。逆時針方向。將其添加到多邊形並繼續，直到返回到初始點。

該算法還將生成一個包含網絡所有外部鍵的大型多邊形。我想你會找到一種方法來識別這個並刪除它。

Answer 2

對於未來的讀者，弗蘭克在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