【问题标题】:Combine nearby polygons合并附近的多边形
【发布时间】:2013-09-08 06:31:45
【问题描述】:

我有一对(闭合的)多边形,每个都定义为一系列点(顶点)。每个多边形代表一块土地,由一条小河隔开,因此溪流在两个多边形之间形成了一个狭窄的间隙。

我正在寻找一种算法来识别和消除间隙,方法是将两个多边形连接成一个连接的多边形。

下图显示了一个示例,其中原始多边形是绿色和红色,生成的多边形显示为黄色。

到目前为止,我已经能够做到以下几点:

  • 对于多边形 A 中的每条边,找到多边形 B 中最近的顶点。
  • 查找多边形 B 在多边形 A 一定距离内的所有顶点。

但我不确定我现在需要做什么。

【问题讨论】:

  • 需要一个更复杂的例子。例如,想象一下将黄色块滑到红色块上,所以你有两个“脚”靠在一个表面上。这样就加入了?或者只有在整个过程中都是一个狭窄的连接?如果它确实被连接起来,你是在中间留下一个整体,还是一个大的连接件?
  • @user1676075:我假设需要对某些阈值进行一些规范。例如,即使在上图中,黄色只是加入它们的一种可能结果。使用更严格的阈值,“北部”“河口”可能太宽,因此那里的连接会更“南部”,导致生成的多边形出现缩进。
  • 你可以缝合顶点

标签: algorithm geometry polygon


【解决方案1】:

使用 Boost.Geometry,您可以调用 buffer 并使用正距离 boost::geometry::strategy::buffer::distance_symmetric<double>{r} 足够大以弥合绿色和红色多边形之间的间隙,然后使用 boost::geometry::strategy::buffer::distance_symmetric<double>{-r} 调用 buffer。如果第二次调用产生伪影,请在r 中添加一个非常小的数字,例如boost::geometry::strategy::buffer::distance_symmetric<double>{-r+0.001}

【讨论】:

  • 如果我们需要自动化这个过程,我们如何找出 r 的值?
【解决方案2】:

为了完整起见,我想分享我的解决方案,作为 python 实现。这是基于 Retsam 接受的答案,以及 David Eisenstat 在 cmets 中对该答案提出的想法,即用来自该多边形的中间顶点替换连接到同一原始多边形上的顶点的凸包上的边.

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)

【讨论】:

    【解决方案3】:

    这是一种粗略的、可扩展的方法。

    1. 将您的空间量化为具有 Z×Z 大小的单元格的网格,并将每个单元格命名为索引 (i,j)。
    2. 对于每个多边形 P,对于每个顶点 V,将 (P,V) 标识为包含它的单元 (i,j)。
    3. 对于每个单元格 (i,j),考虑一组多边形-顶点对 (P_k,V_k),k = 1...K,已用它识别。合并多边形 P_a 和 P_b 的顶点 V_a 和 V_b 当且仅当 P_a 和 P_b 不是同一个多边形并且 V_a 和 V_b 是这些对中最接近的。重复合并顶点,直到没有更多要合并的顶点。合并的顶点获取源顶点位置的平均值。

    【讨论】:

      【解决方案4】:

      你可以试试morphological operations。具体来说,您可以尝试膨胀,然后是侵蚀(也称为形态“闭合”)。 n 个像素的膨胀(其中 n 大于河流的宽度)将组合形状。随后的侵蚀将消除对图形其余部分造成的大部分损害。它不会是完美的(它完美地结合这两个形状,但代价是形状的其余部分的一些软化),但也许从整体操作的结果来看,你可以想想办法改正它。

      一般来说,这些形态学操作是在位图上完成的,而不是在多边形上。但是在多边形的角上运行简单的操作可能会起作用。

      【讨论】:

      • 这是一个有趣的想法,但我认为这对我的目的来说太有损了。但总是乐于了解新事物。
      【解决方案5】:

      您可能会查看凸包算法的修改。凸包算法采用一组点并绘制包含这些点的最小凸形状。您的问题几乎是一个凸包问题,除了顶部的那些凹入区域。简单地使用凸包算法会给你这个,这是接近的,但不是你所需要的(注意。不同的棕色区域)

      根据您要执行的操作,凸包可能“足够好”,但如果不是,您仍然可以修改算法以忽略非凸部分,然后简单地合并两个多边形。

      特别是 this pdf 显示了如何合并两个凸包,这很像您正在尝试做的事情。

      【讨论】:

      • 感谢您的想法。我正在考虑使用凸包,但我正在使用的一些形状非常凹,以至于凸包会被淘汰。但也许我可以像你建议的那样使用 CH 算法作为起点。
      • 是的;具体来说,我想说看看“合并”或“分而治之”算法,因为它专门处理组合现有船体。]
      • @sh1ftst0rm 如果您采用凸包并将连接来自同一多边形的顶点的边替换为原始边的序列,您会得到黄色区域。我不确定这种方法的所有极端情况是什么。
      • @DavidEisenstat:这是一个非常有趣的想法,我会研究一下。谢谢!
      • @DavidEisenstat:您替换连接同一多边形顶点的 CH 边的想法非常有效。如果您想将其发布为答案,我很乐意投票。
      猜你喜欢
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2017-10-16
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2012-12-28
      相关资源
      最近更新 更多