分割中间多边形之间的交点

Question

我正在处理一些生物成像样本，并尝试创建细胞形状的数字模型。为了简单起见，我想通过将它们建模为多边形来概括它们的形状。

我正在努力将两个重叠的多边形分割成在其交叉点处不共享重叠区域的多边形。相反，该区域分为两个形状。下面最能说明我的意图。

我使用 Python 和 OpenCV 包工作，但很乐意实现任何可以解决此问题的替代包。 （numpythonic 方式是最好的 - 如果可能的话！）

Answer 1

from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
import cv2


threshold = 200
thickness = 16
linewidth = 2
white = 255
black = 0


def fun(image):
    res = np.full(image.shape, white).astype(np.uint8)
    # determine contours
    grayscale = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    thresh = cv2.threshold(grayscale, threshold, white, cv2.THRESH_BINARY)[1]
    contours = cv2.findContours(thresh, mode=cv2.RETR_TREE, method=cv2.CHAIN_APPROX_NONE)[0]
    # create intersection and union
    intersection, union = sorted(contours, key=len)[::len(contours) - 2]
    outline = np.full(res.shape, white).astype(np.uint8)
    cv2.drawContours(outline, [union], -1, black, thickness, cv2.LINE_AA)
    inline = np.full(res.shape, white).astype(np.uint8)
    cv2.drawContours(inline, [intersection], -1, black, thickness, cv2.LINE_AA)
    # determine points for line
    points = np.logical_and(outline == 0, inline == 0).astype(np.uint8)[:,:,0] * white
    a, b = cv2.findContours(points, mode=cv2.RETR_TREE, method=cv2.CHAIN_APPROX_NONE)[0]
    a = np.mean(a.squeeze(), axis=0).astype(int)
    b = np.mean(b.squeeze(), axis=0).astype(int)
    # draw outline of union
    cv2.drawContours(res, [union], -1, black, linewidth, cv2.LINE_AA)
    # draw connecting line
    cv2.line(res, a, b, black, linewidth)
    return res

for img in ['img1.jpg', 'img2.jpg']:
    orig = np.asarray(Image.open(img))
    res = fun(orig)
    fig, ax = plt.subplots(1, 2)
    for i, img in enumerate([orig, res]):
        ax[i].imshow(img, cmap='gray');
        ax[i].axes.get_xaxis().set_visible(False)
        ax[i].axes.get_yaxis().set_visible(False)
    plt.show();

此方法的工作原理如下：

1. 首先将图像转换为灰度并使用阈值对图像进行二值化，然后确定

   contours

   。接下来，通过增加大小对轮廓进行排序，这确保交集的轮廓（要删除的部分）排在第一位，并集的轮廓（要保留的部分）排在倒数第二位（最后是图像边界的轮廓）。
2. 创建两张新的白色图像，在一张上用黑色绘制并集

   outline

   的轮廓，在另一张上用黑色绘制交集

   inline

   的轮廓，然后将两个图像的黑色部分相交以获得标记端点的补丁对于新线路。然后计算新线的端点

   a

   和

   b

   作为它们各自补丁的平均值。
3. 在新的、最初为白色的黑色图像上绘制一条连接两点

   a

   和

   b

   的直线以及联合的轮廓，从而得到所需的输出图像。

备注：这种方法不适用于第三种情况（因为它假设形状是凸的，或者更确切地说，只有一个交叉点），但我的猜测是，对其进行调整以适用于这种情况应该是可行的。

Answer 2

shapely是一个提供各种几何操作的库。

以下代码示例根本不适合生产，仍然需要相当多的微调才能涵盖所有可能的情况（例如，星形的结果还不理想），但它说明了形状功能如何帮助构建这样的功能：

import shapely
import shapely.ops
from matplotlib import pyplot as plt
from shapely.geometry import box, LineString, Polygon, Point
from shapely import plotting as plotter

def distribute_intersections(poly1, poly2, grid_size=0.0001):
    intersection_points = poly1.boundary.intersection(poly2.boundary, grid_size=grid_size)
    intersection_line = LineString(shapely.get_coordinates(intersection_points))
    intersections = poly1.intersection(poly2, grid_size=grid_size)
    split_lines = shapely.intersection(intersections, intersection_line, grid_size=grid_size)

    # Loop through all split lines found and apply them one by one
    poly1_distributed = poly1
    poly2_distributed = poly2
    for split_line in shapely.get_parts(split_lines):
        if isinstance(split_line, Point):
            continue

        # For both polygons, split them with the current split line and remove the part
        # that is ~covered by the other polygon
        poly1_distributed = shapely.ops.split(poly1_distributed, split_line)
        poly1_distributed = _remove_covered_parts(poly1_distributed, poly2, grid_size)
        poly2_distributed = shapely.ops.split(poly2_distributed, split_line)
        poly2_distributed = _remove_covered_parts(poly2_distributed, poly1, grid_size)

    return poly1_distributed, poly2_distributed

def _remove_covered_parts(poly_split, intersecting_poly, grid_size):
    for poly in shapely.get_parts(poly_split):
        intersection = poly.intersection(intersecting_poly, grid_size=grid_size)
        if intersection.area > 0.9 * poly.area:
            continue

        return poly

poly1 = box(0, 0, 10, 10)
poly2 = box(7, 7, 15, 20)
poly1_result, poly2_result = distribute_intersections(poly1, poly2)

circle1 = Point(0, 5).buffer(6, resolution=4)
circle2 = Point(0, 12).buffer(6, resolution=4)
circle1_result, circle2_result = distribute_intersections(circle1, circle2)

circle3 = Point(5, 0).buffer(6, resolution=4)
star1 = Polygon([
    (8, 0), (12, -1), (13, -5), (14, -1), (20, 0), (16, 1), (18, 6), (14, 4),
    (8, 3), (12, 2), (8, 0)
])
circle3_result, star1_result = distribute_intersections(circle3, star1)

# Plot results
_, ax = plt.subplots(ncols=3)
plotter.plot_polygon(poly1_result, ax=ax[0], color="red", alpha=0.1)
plotter.plot_polygon(poly2_result, ax=ax[0], color="green", alpha=0.1)

plotter.plot_polygon(circle1_result, ax=ax[1], color="red", alpha=0.1)
plotter.plot_polygon(circle2_result, ax=ax[1], color="green", alpha=0.1)

plotter.plot_polygon(circle3_result, ax=ax[2], color="red", alpha=0.1)
plotter.plot_polygon(star1_result, ax=ax[2], color="green", alpha=0.1)

for cur_ax in ax:
    cur_ax.set_aspect("equal")
plt.show()

结果：