2013-09-04 3 views
3

У меня есть пара (замкнутых) многоугольников, каждая из которых определена как последовательность точек (вершин). Каждый из полигонов представляет собой участок земли, разделенный небольшой рекой, поэтому поток образует узкую щель между двумя полигонами.Объединить близлежащие полигоны

Я ищу алгоритм для выявления и устранения пробела путем объединения двух полигонов в один многоугольник.

На рисунке ниже показан пример, где исходные многоугольники являются зеленым и красным, а полученный многоугольник показан желтым цветом.

example

До сих пор я был в состоянии сделать следующее:

  • Для каждого ребра многоугольника-A, найти ближайшую вершину многоугольника-B.
  • Найти все вершины многоугольника-B, которые находятся на определенном расстоянии от полигона-A.

Но я не совсем уверен, что мне нужно делать сейчас.

+0

Нужен более сложный пример. Например, представьте, что скользящий желтый блок против красного цвета, так что у вас есть две «ноги» на одной поверхности. Соединяется ли это? Или только если это узкая связь в целом? Если он соединяется, вы остаетесь с целым посередине или с одной большой частью? – user1676075

+0

@ user1676075: Я предполагаю, что потребуется некоторая спецификация некоторых пороговых значений. Например, даже на рисунке выше желтый является всего лишь одним из возможных результатов их присоединения. С более строгими порогами «северный» «устье реки» может быть слишком широким, и поэтому соединение там будет дальше «на юг», что приведет к отступу в полученном многоугольнике. – brianmearns

+0

Вы можете шить вершины –

ответ

3

Вы можете посмотреть модификацию алгоритма выпуклого корпуса. Алгоритм выпуклого корпуса принимает множество точек и рисует минимальную выпуклую форму, содержащую эти точки. Ваш вопрос - это почти выпуклый вопрос корпуса, за исключением тех вогнутых областей наверху. Просто с помощью алгоритма оконтуривающего даст вам это, что близко, но не совсем то, что вам нужно (обратите внимание .the коричневых участков, которые отличаются)

Convex Hull

В зависимости от того, что вы пытаетесь сделать, выпуклый корпус может быть «достаточно хорош», но если нет, вы все равно сможете модифицировать алгоритм, чтобы игнорировать части, которые не являются выпуклыми, а simpy объединяет два многоугольника.

В частности, this pdf показывает, как можно объединить два выпуклых корпуса, что очень похоже на то, что вы пытаетесь сделать.

+0

Спасибо за идею. Я рассматривал возможность использования выпуклой оболочки, но некоторые из форм, с которыми я работаю, настолько вогнуты, что выпуклый корпус будет уходить. Но, возможно, я могу использовать алгоритм CH как отправную точку, как вы предлагаете. – brianmearns

+0

Да; в частности, я бы сказал, заглядывая в алгоритм «слияния» или «делить и побеждать», поскольку он конкретно связан с объединением существующих корпусов.] – Retsam

+2

@ sh1ftst0rm Вы получили бы желтую область, если бы вы взяли выпуклый корпус и заменили ребра соединяющий вершины из того же многоугольника с последовательностью ребер в оригинале. Я не уверен, что все из этих угловых случаев этого подхода. –

1

Вы можете попробовать morphological operations. В частности, вы можете попробовать дилатацию, а затем erosion (также известный как морфологический «закрытие»). Расширение на n пикселей - где n больше, чем ширина реки - объединит формы. Последующая эрозия затем отменит большую часть ущерба, нанесенного остальной части фигуры. Это было бы не идеально (это было бы прекрасно сочетать две формы, но за счет некоторого смягчения остальной формы), но, возможно, с результатом полной операции, вы могли бы выяснить способ исправить Это.

Обычно эти морфологические операции выполняются на растровых изображениях, а не на многоугольниках. Но выполнение простых операций над углами полигонов может работать.

+0

Это интересная идея, но я думаю, что это будет слишком скучно для моих целей. Но всегда рады узнать о чем-то новом. – brianmearns

0

Это приблизительный, масштабируемый подход.

  1. Пронумеруйте свое пространство в сетке с ячейками размера Z-by-Z и назовите каждую ячейку своими индексами (i, j).
  2. Для каждого многоугольника P для каждой вершины V идентифицируйте (P, V) с содержащейся в нем ячейкой (i, j).
  3. Для каждой ячейки (i, j) рассмотрим множество многоугольных вершинных пар (P_k, V_k), k = 1 ... K, которые были идентифицированы с ним. Сверните вершины V_a и V_b многоугольников P_a и P_b тогда и только тогда, когда P_a и P_b не являются одним и тем же многоугольником, а V_a и V_b являются ближайшими среди таких пар. Повторяйте слияние вершин до тех пор, пока больше не будет слияния. Слитые вершины получают среднее положение позиций вершин.
5

Для полноты, я хотел поделиться своим решением, как реализация python. Это основано на принятом ответе Retsam и идее, представленной Дэвидом Эйзенстатом в комментариях к этому ответу, который заключается в замене ребер на выпуклой оболочке, которые соединяются с вершинами на одном и том же исходном многоугольнике, с промежуточными вершинами из этого многоугольника ,

def joinPolygons(polya, polyb): 
    """ 
    Generate and return a single connected polygon which includes the two given 
    polygons. The connection between the two polygons is based on the convex hull 
    of the composite polygon. All polygons are sequences of two-tuples giving the 
    vertices of the polygon as (x, y), in order. That means vertices that are adjacent 
    in the sequence are adjacent in the polygon (connected by an edge). The first and 
    last vertices in the sequence are also connected by any edge (implicitly closed, do 
    not duplicate the first point at the end of the sequence to close it). 

    Only simple polygons are supported (no self-intersection). 
    """ 

    #Just to make it easier to identify and access by index. 
    polygons = [polya, polyb] 

    #Create a single list of points to create the convex hull for (each 
    # point is a vertex of one of the polygons). 
    #Additionally, each point includes some additional "cargo", indicating which 
    # polygon it's from, and it's index into that polygon 
    # This assumes the implementation of convexHull simply ignores everything 
    # beyond the first two elements of each vertex. 
    composite = [] 
    for i in range(len(polygons)): 
     points = polygons[i] 
     composite += [(points[j][0], points[j][1], j, i) for j in xrange(len(points))] 

    #Get the convex hull of the two polygons together. 
    ch = convexHull(composite) 

    #Now we're going to walk along the convex hull and find edges that connect two vertices 
    # from the same source polygon. We then replace that edge with all the intervening edges 
    # from that source polygon. 

    #Start with the first vertex in the CH. 
    x, y, last_vnum, last_pnum = ch[0] 

    #Here is where we will collect the vertices for our resulting polygon, starting with the 
    # first vertex on the CH (all the vertices on the CH will end up in the result, plus some 
    # additional vertices from the original polygons). 
    results = [(x, y)] 

    #The vertices of the convex hull will always walk in a particular direction around each 
    # polygon (i.e., forwards in the sequence of vertices, or backwards). We will use this 
    # to keep track of which way they go. 
    directions = [None for poly in polygons] 

    #Iterate over all the remaining points in the CH, and then back to the first point to 
    # close it. 
    for x, y, vnum, pnum in list(ch[1:]) + [ch[0]]: 

     #If this vertex came from the same original polygon as the last one, we need to 
     # replace the edge between them with all the intervening edges from that polygon. 
     if pnum == last_pnum: 

      #Number of vertices in the polygon 
      vcount = len(polygons[pnum]) 

      #If an edge of the convex hull connects the first and last vertex of the polygon, 
      # then the CH edge must also be an edge of the polygon, because the two vertices are 
      # adjacent in the original polygon. Therefore, if the convex 
      # hull goes from the first vertex to the last in a single edge, then it's walking 
      # backwards around the polygon. Likewise, if it goes from the last to the first in 
      # a single edge, it's walking forwards. 
      if directions[pnum] is None: 
       if last_vnum < vnum: 
        if last_vnum == 0 and vnum == vcount - 1: 
         direction = -1 
        else: 
         direction = 1 
       else: 
        if last_vnum == vcount - 1 and vnum == 0: 
         direction = 1 
        else: 
         direction = -1 
       directions[pnum] = direction 
      else: 
       direction = directions[pnum] 

      #Now walk from the previous vertex to the current one on the source 
      # polygon, and add all the intevening vertices (as well as the current one 
      # from the CH) onto the result. 
      v = last_vnum 
      while v != vnum: 
       v += direction 
       if v >= vcount: 
        v = 0 
       elif v == -1: 
        v = vcount - 1 
       results.append(polygons[pnum][v]) 

     #This vertex on the CH is from a different polygon originally than the previous 
     # vertex, so we just leave them connected. 
     else: 
      results.append((x, y)) 

     #Remember this vertex for next time. 
     last_vnum = vnum 
     last_pnum = pnum 

    return results 



def convexHull(points, leftMostVert=None): 
    """ 
    Returns a new polygon which is the convex hull of the given polygon. 

    :param: leftMostVert The index into points of the left most vertex in the polygon. 
          If you don't know what it is, pass None and we will figure it 
          out ourselves. 
    """ 
    point_count = len(points) 

    #This is implemented using the simple Jarvis march "gift wrapping" algorithm. 
    # Generically, to find the next point on the convex hull, we find the point 
    # which has the smallest clockwise-angle from the previous edge, around the 
    # last point. We start with the left-most point and a virtual vertical edge 
    # leading to it. 

    #If the left-most vertex wasn't specified, find it ourselves. 
    if leftMostVert is None: 
     minx = points[0][0] 
     leftMostVert = 0 
     for i in xrange(1, point_count): 
      x = points[i][0] 
      if x < minx: 
       minx = x 
       leftMostVert = i 

    #This is where we will build up the vertices we want to include in the hull. 
    # They are stored as indices into the sequence `points`. 
    sel_verts = [leftMostVert] 

    #This is information we need about the "last point" and "last edge" in order to find 
    # the next point. We start with the left-most point and a pretend vertical edge. 

    #The index into `points` of the last point. 
    sidx = leftMostVert 

    #The actual coordinates (x,y) of the last point. 
    spt = points[sidx] 

    #The vector of the previous edge. 
    # Vectors are joined tail to tail to measure angle, so it 
    # starts at the last point and points towards the previous point. 
    last_vect = (0, -1, 0) 
    last_mag = 1.0 

    #Constant 
    twopi = 2.0*math.pi 

    #A range object to iterate over the vertex numbers. 
    vert_nums = range(point_count) 

    #A list of indices of points which have been determined to be colinear with 
    # another point and a selected vertex on the CH, and which are not endpoints 
    # of the line segment. These points are necessarily not vertices of the convex 
    # hull: at best they are internal to one of its edges. 
    colinear = [] 

    #Keep going till we come back around to the first (left-most) point. 
    while True: 
     #Try all other end points, find the one with the smallest CW angle. 
     min_angle = None 
     for i in vert_nums: 

      #Skip the following points: 
      # -The last vertex (sidx) 
      # -The second to last vertex (sel_verts[-2]), that would just backtrack along 
      # the edge we just created. 
      # -Any points which are determined to be colinear and internal (indices in `colinear`). 
      if i == sidx or (len(sel_verts) > 1 and i == sel_verts[-2]) or i in colinear: 
       continue 

      #The point to test (x,y) 
      pt = points[i] 

      #vector from current point to test point. 
      vect = (pt[0] - spt[0], pt[1] - spt[1], 0) 
      mag = math.sqrt(vect[0]*vect[0] + vect[1]*vect[1]) 

      #Now find clockwise angle between the two vectors. Start by 
      # finding the smallest angle between them, using the dot product. 
      # Then use cross product and right-hand rule to determine if that 
      # angle is clockwise or counter-clockwise, and adjust accordingly. 

      #dot product of the two vectors. 
      dp = last_vect[0]*vect[0] + last_vect[1]*vect[1] 
      cos_theta = dp/(last_mag * mag) 

      #Ensure fp erros don't become domain errors. 
      if cos_theta > 1.0: 
       cos_theta = 1.0 
      elif cos_theta < -1.0: 
       cos_theta = -1.0 

      #Smaller of the two angles between them. 
      theta = math.acos(cos_theta) 

      #Take cross product of last vector by test vector. 
      # Except we know that Z components in both input vectors are 0, 
      # So the X and Y components of the resulting vector will be 0. Plus, 
      # we only care aboue the Z component of the result. 
      cpz = last_vect[0]*vect[1] - last_vect[1]*vect[0] 

      #Assume initially that angle between the vectors is clock-wise. 
      cwangle = theta 
      #If the cross product points up out of the plane (positive Z), 
      # then the angle is actually counter-clockwise. 
      if cpz > 0: 
       cwangle = twopi - theta 

      #If this point has a smaller angle than the others we've considered, 
      # choose it as the new candidate. 
      if min_angle is None or cwangle < min_angle: 
       min_angle = cwangle 
       next_vert = i 
       next_mvect = vect 
       next_mag = mag 
       next_pt = pt 

      #If the angles are the same, then they are colinear with the last vertex. We want 
      # to pick the one which is furthest from the vertex, and put all other colinear points 
      # into the list so we can skip them in the future (this isn't just an optimization, it 
      # appears to be necessary, otherwise we will pick one of the other colinear points as 
      # the next vertex, which is incorrect). 
      #Note: This is fine even if this turns out to be the next edge of the CH (i.e., we find 
      # a point with a smaller angle): any point with is internal-colinear will not be a vertex 
      # of the CH. 
      elif cwangle == min_angle: 
       if mag > next_mag: 
        #This one is further from the last vertex, so keep it as the candidate, and put the 
        # other colinear point in the list. 
        colinear.append(next_vert) 
        min_angle = cwangle 
        next_vert = i 
        next_mvect = vect 
        next_mag = mag 
        next_pt = pt 
       else: 
        #This one is closer to the last vertex than the current candidate, so just keep that 
        # as the candidate, and put this in the list. 
        colinear.append(i) 

     #We've found the next vertex on the CH. 
     # If it's the first vertex again, then we're done. 
     if next_vert == leftMostVert: 
      break 
     else: 
      #Otherwise, add it to the list of vertices, and mark it as the 
      # last vertex. 
      sel_verts.append(next_vert) 
      sidx = next_vert 
      spt = next_pt 
      last_vect = (-next_mvect[0], -next_mvect[1]) 
      last_mag = next_mag 

    #Now we have a list of vertices into points, but we really want a list of points, so 
    # create that and return it. 
    return tuple(points[i] for i in sel_verts) 
+0

это потрясающе! Вы просто спасли мне кучу работы! –

Смежные вопросы