python - 如何将匀称的几何图形投影到不同的形状

Question

我有一个 GeoPandas GeoDataFrame 形式的模板，其中包含与页面上的区域相对应的几个形状优美的几何图形。

对于这个过程，页面被拍照，我需要模板中的几何图形以及拍照页面上的相应区域。现在，我将模板“投影”到拍摄页面的边界框与拍摄的图像不对齐。我也很确定有更好的方法来做到这一点。

我创建了一个简化的示例来说明：

我有这个文档，两个轴都缩放为 1。

我创建了与图像上的对象匹配的模板

shapes = [
    ('line1', LineString([(0.146, 0.216), (0.476 , 0.216)])),
    ('line2', LineString([(0.498, 0.871), (0.838, 0.871)]))
]

shapes = gpd.GeoDataFrame(shapes, columns=['name', 'geometry'], geometry='geometry')

将其绘制在图像顶​​部，我们可以验证形状是否与图像中的形状紧密对齐。

plt.imshow(img, extent=[0,1,1,0], aspect=1.4142)
ax = plt.gca()
shapes.plot(ax=ax)
ax.set_aspect(1.4142)

我为同一个文档拍了一张照片，得到了页面的边界框。

photo = load_image('example_photo.jpg')

bbox = Polygon([(0.847429096698761, 0.047594085335731506),
                (0.9442346692085266, 0.8787651062011719),
                (0.05563090369105339, 0.8789013028144836),
                (0.12740036845207214, 0.06380380690097809),
                (0.847429096698761, 0.047594085335731506)])

plt.figure(figsize=(6,6))
plt.imshow(photo, extent=[0,1,1,0], aspect=1.4142)
plt.plot(*bbox.boundary.xy)

问题出在下一步，因为我尝试将原始模板重新映射或投影到边界框的形状中。这是我尝试过的，但我确信这不是最有效的方法。它也不起作用。

下面的方法总结如下：

1. 找出哪个边缘是顶部，底部，右侧，左侧并定位它们。
2. 将每个点映射到新形状中，找到与边界框的每个相对侧相交的线的交点，与单位正方形中的相同位置成比例。

def get_edges(bbox, visualize=True):
    '''Takes shapely polygon with exactly 5 points'''

    # check for 5 points
    if len(bbox.boundary.coords) != 5:
        raise('Polygon must have 5 points (4 sides)')
            
    #find top/bottom edge
    x, y = bbox.boundary.xy
    
    # remove last point which must be the same as the first
    x = np.array(x)[:-1]
    y = np.array(y)[:-1]
    
    # sort by y values.  Higher is closer to top, lower closer to bottom
    y_sorted = np.argsort(y)
    
    # get the index of the top and bottom lines
    top_points_idx = y_sorted[-2:]
    bot_points_idx = y_sorted[:2]
    
    # order the top point coords left to right
    top_point_order = top_points_idx[
        np.argsort(x[top_points_idx])
    ]
    
    bot_point_order = bot_points_idx[
        np.argsort(x[bot_points_idx])
    ]
        
    
    top_points = np.array(bbox.boundary.coords)[top_point_order]
    bot_points = np.array(bbox.boundary.coords)[bot_point_order]
    
    left_points = LineString([bot_points[0], top_points[0]])
    right_points = LineString([bot_points[1], top_points[1]])
    
    top_points = LineString(top_points)
    bot_points = LineString(bot_points)
    
    return top_points, bot_points, left_points, right_points

def project_unit_square_point(point, bbox, visualize=False):

    # check for 5 points
    if len(bbox.boundary.coords) != 5:
        raise('Polygon must have 5 points (4 sides)')

    # use the position as the portion of the each side
    x_scale, y_scale = point.coords[0]
    
    top, bot, left, right = get_edges(bbox)

    # find proportional intersections on edges
    top_point = top.interpolate(x_scale*top.length)
    bot_point = bot.interpolate(x_scale*bot.length)
    left_point = left.interpolate(y_scale*left.length)
    right_point = right.interpolate(y_scale*right.length)
    
    # connect edge points
    vline = LineString([top_point, bot_point])
    hline = LineString([left_point, right_point])
    
    # new point is intersection of vline and hline
    new_point = vline.intersection(hline)
    
    return new_point

def project_unit_square_geom(geom, bbox):
    
    new_points = []
    for point in geom.coords:
        new_points.append(project_unit_square_point(Point(point), bbox))
        
    new_geom = LineString(new_points)
    
    return new_geom

# project geoms onto form
projected_shapes = []
for shape in shapes.geometry:
    projected_shapes.append(
        project_unit_square_geom(shape, bbox)
    )

# create a new df for the mapped shapes
projected_shapes = gpd.GeoSeries(projected_shapes, name='geometry')

projected_shapes = gpd.GeoDataFrame({'name': shapes['name'],
                                     'geometry': projected_shapes},
                                    geometry='geometry')

然后，当我可视化结果时，我得到了这个：

plt.figure(figsize=(6,6))
plt.imshow(photo, extent=[0,1,1,0])
plt.plot(*bbox.boundary.xy)
ax = plt.gca()
projected_shapes.plot(ax=ax)
ax.set_aspect(1.4142)

接近但不够接近。Obvioulsy我的方法不起作用。如何将模板形状映射到边界框定义的新形状上？

这是要使用的原始图像。

标签： pythongeopandasshapely