我正在处理一些生物成像样本,并尝试创建细胞形状的数字模型。为了简单起见,我想通过将它们建模为多边形来概括它们的形状。
我正在努力将两个重叠的多边形分割成在其交叉点处不共享重叠区域的多边形。相反,该区域分为两个形状。下面最能说明我的意图。
我使用 Python 和 OpenCV 包工作,但很乐意实现任何可以解决此问题的替代包。 (numpythonic 方式是最好的 - 如果可能的话!)
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();
此方法的工作原理如下:
contours
。接下来,通过增加大小对轮廓进行排序,这确保交集的轮廓(要删除的部分)排在第一位,并集的轮廓(要保留的部分)排在倒数第二位(最后是图像边界的轮廓)。 outline
的轮廓,在另一张上用黑色绘制交集inline
的轮廓,然后将两个图像的黑色部分相交以获得标记端点的补丁对于新线路。然后计算新线的端点 a
和 b
作为它们各自补丁的平均值。a
和 b
的直线以及联合的轮廓,从而得到所需的输出图像。备注:这种方法不适用于第三种情况(因为它假设形状是凸的,或者更确切地说,只有一个交叉点),但我的猜测是,对其进行调整以适用于这种情况应该是可行的。
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()
结果: