Полигоны из сети подключенных точек

Question

Учитывая массив двумерных точек (#pts x 2) и массив из которых связаны с которыми (#bonds x 2 int array с индексами 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