找到不规则形状多边形的"视觉"中心的最快方法是什么？

Question

我需要找到一个点,它是一个不规则形状的多边形的视觉中心.通过视觉中心,我的意思是在视觉上看起来位于多边形的大区域的中心的点.应用程序是在多边形内部放置一个标签.

这是一个使用内部缓冲的解决方案:

https://web.archive.org/web/20150708063910/http://proceedings.esri.com/library/userconf/proc01/professional/papers/pap388/p388.htm

如果要使用它,找到缓冲区的有效且快速的方法是什么？如果要使用任何其他方式,这是哪种方式？

真正坚韧的多边形的一个很好的例子是一个巨大的厚U(用Arial Black或Impact或一些这样的字体书写).

Answer 1

我从MapBox中找到了一个非常好的解决方案,名为Polylabel.完整的源代码也可以在他们的Github上找到.

基本上它试图找到多边形的视觉中心,如T奥斯汀所说.

某些细节表明这可能是一个实用的解决方案:

不幸的是,计算[理想解决方案]既复杂又缓慢.已发布的问题解决方案需要约束Delaunay三角剖分或计算直线骨架作为预处理步骤 - 这两个步骤都很慢且容易出错.

对于我们的用例,我们不需要一个确切的解决方案 - 我们愿意交易一些精度以获得更快的速度.当我们在地图上放置标签时,以毫秒为单位计算它比在数学上完美更重要.

关于使用的快速说明.源代码适用于开箱即用的Javascript,但是如果你打算使用"普通"多边形,那么你应该将它包装在一个空数组中,因为这里的函数采用GeoJSONPolygons而不是普通多边形,即

var myPolygon = [[x1, y1], [x2, y2], [x3, y3]];
var center = polylabel([myPolygon]);

Answer 2

如果可以将多边形转换为二进制图像,则可以使用图像处理领域中存在的基础,例如:块代表二进制图像的快速骨架算法.

但由于离散化错误和额外工作,这在一般情况下并不合理.

但是,也许您会发现这些有用:

• 简单多边形的直骨架
• 在(几乎)线性时间中确定简单多边形的骨架

编辑:也许你想要寻找多边形中包含的最大圆的中心点.它不一定总是在观察中心,但大多数时间可能会给出预期的结果,并且只在略微病态的情况下才能完全关闭.

Answer 3

怎么样:

如果多边形的质心在多边形内,则使用它,否则:

1)从质心通过多边形延伸一条线,将多边形分成两半相等的区域

2)"视觉中心"是在线接触周边的最近点和在远离质心的方向上切割周边的下一点之间的点的中间点

这里有几张图片来说明它:

Answer 4

您是否考虑过使用质心公式？

http://en.wikipedia.org/wiki/Centroid

http://en.wikipedia.org/wiki/K-means_algorithm

Answer 5

质心法已经被多次建议。我认为这是一个很好的资源，它非常直观地描述了该过程（以及许多其他有用的多边形技巧）：

http://paulbourke.net/geometry/polygonmesh/centroid.pdf

另外，为了放置一个简单的 UI 标签，只需计算多边形的边界框（由多边形中任何顶点的最低和最高 x 和 y 坐标定义的矩形）并使其中心位于：

{
    x = min_x + (max_x - min_x)/2,
    y = min_y + (max_y - min_y)/2
}

这比计算质心要快一些，这对于实时或嵌入式应用程序可能很重要。

另请注意，如果您的多边形是静态的（它们不会改变形式），您可以通过将 BB 中心/质心计算的结果（相对于例如多边形的第一个顶点）保存到以下数据结构来进行优化多边形。

Answer 6

以下是我尝试过的六种不同方法。

\n

\n
1. cv2基于质心 ( get_center_of_mass)
\n
2. shapely基代表点 ( get_representative_point)
\n
3. cv2+基于骨架skimage.skeleton形状的质心( )get_skeleton_center_of_mass
\n
4. scipy基于到边界的最远距离 ( get_furthest_point_from_edge)
\n
5. cv2基于先前最远边界距离算法的版本 ( get_furthest_point_from_edge_cv2)
\n
6. get_center_of_half_area_line@T.Austin ( )在该线程中提出的“半区域线的中心点”算法
\n

\n

让我们从导入和一些辅助函数开始

\n

import numpy as np\nimport cv2\nfrom shapely.geometry import Polygon, LineString, MultiLineString, Point, MultiPoint, GeometryCollection\nfrom skimage.morphology import skeletonize, medial_axis\nfrom scipy.optimize import minimize_scalar\nfrom scipy.ndimage.morphology import distance_transform_edt\nimport matplotlib.pyplot as plt\n\nH, W = 300, 300\n\ndef get_random_contour():\n    xs = np.random.randint(0, W, 4)\n    ys = np.random.randint(0, H, 4)\n    cnt = np.array([[x,y] for x,y in zip(xs,ys)])\n    mask = draw_contour_on_mask((H,W), cnt)\n    cnt, _ = cv2.findContours(mask, 1, 2)\n    cnt = cnt[0]\n    return cnt\n\ndef draw_contour_on_mask(size, cnt, color:int = 255):\n    mask = np.zeros(size, dtype=\'uint8\')\n    mask = cv2.drawContours(mask, [cnt], -1, color, -1)\n    return mask\n\ndef get_center_of_mass(cnt):\n    M = cv2.moments(cnt)\n    cx = int(M[\'m10\']/M[\'m00\'])\n    cy = int(M[\'m01\']/M[\'m00\'])\n    return cx, cy\n\ndef split_mask_by_line(mask, centroid:tuple, theta_degrees:float, eps:float = 1e-4):\n    h, w = mask.shape[:2]\n    mask = mask[..., None]\n    cx, cy = centroid\n    # convert theta first to radians and then to line slope(s)\n    theta_degrees = np.atleast_1d(theta_degrees)\n    theta_degrees = np.clip(theta_degrees, -90+eps, 90-eps)\n    theta_rads = np.radians(theta_degrees)\n    slopes = np.tan(theta_rads)[:, None]\n    # define the line(s)\n    x = np.arange(w, dtype="int32")\n    y = np.int32(slopes * (x - cx) + cy)\n    _y = np.arange(h, dtype="int32")\n    # split the input mask into two halves by line(s)\n    m = (y[..., None] <= _y).T\n    m1 = (m * mask).sum((0,1))\n    m2 = ((1 - m) * mask).sum((0,1))\n    m2 = m2 + eps if m2==0 else m2\n    # calculate the resultant masks ratio\n    ratio = m1/m2\n    return (x.squeeze(), y.squeeze()), ratio\n\ndef get_half_area_line(mask, centroid: tuple, eps: float = 1e-4):\n    # find the line that splits the input mask into two equal area halves\n    minimize_fun = lambda theta: abs(1. - split_mask_by_line(mask, centroid, theta, eps=eps)[1].item())\n    bounds = np.clip((-90, 90), -90 + eps, 90 - eps)\n    res = minimize_scalar(minimize_fun, bounds=bounds, method=\'bounded\')\n    theta_min = res.x\n    line, _ = split_mask_by_line(mask, centroid, theta_min)\n    return line\n

\n

现在让我们定义寻找视觉中心的函数

\n

def get_representative_point(cnt):\n    poly = Polygon(cnt.squeeze())\n    cx = poly.representative_point().x\n    cy = poly.representative_point().y\n    return cx, cy\n\ndef get_skeleton_center_of_mass(cnt):\n    mask = draw_contour_on_mask((H,W), cnt)\n    skel = medial_axis(mask//255).astype(np.uint8) #<- medial_axis wants binary masks with value 0 and 1\n    skel_cnt,_ = cv2.findContours(skel,1,2)\n    skel_cnt = skel_cnt[0]\n    M = cv2.moments(skel_cnt) \n    if(M["m00"]==0): # this is a line\n        cx = int(np.mean(skel_cnt[...,0]))\n        cy = int(np.mean(skel_cnt[...,1]))\n    else:\n        cx = int(M[\'m10\']/M[\'m00\'])\n        cy = int(M[\'m01\']/M[\'m00\'])\n    return cx, cy\n\ndef get_furthest_point_from_edge(cnt):\n    mask = draw_contour_on_mask((H,W), cnt)\n    d = distance_transform_edt(mask)\n    cy, cx = np.unravel_index(d.argmax(), d.shape)\n    return cx, cy\n\ndef get_furthest_point_from_edge_cv2(cnt):\n    mask = draw_contour_on_mask((H,W), cnt)\n    dist_img = cv2.distanceTransform(mask, distanceType=cv2.DIST_L2, maskSize=5).astype(np.float32)\n    cy, cx = np.where(dist_img==dist_img.max())\n    cx, cy = cx.mean(), cy.mean()  # there are sometimes cases where there are multiple values returned for the visual center\n    return cx, cy\n\ndef get_center_of_half_area_line(cnt):\n    mask = draw_contour_on_mask((H,W), cnt, color=1)\n    # get half-area line that passes through centroid\n    cx, cy = get_center_of_mass(mask)\n    line = get_half_area_line(mask, centroid=(cx, cy))\n    line = LineString(np.array(list(zip(line))).T.reshape(-1, 2))\n    # find the visual center\n    contours, _ = cv2.findContours(mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)\n    contours = [c for c in contours if cv2.contourArea(c) > 5]\n    polys = [Polygon(c.squeeze(1)) for c in contours if len(c) >= 3]  # `Polygon` must have at least 3 points\n    cpoint = Point(cx, cy)\n    points = []\n    for poly in polys:\n        isect = poly.intersection(line)\n        if isect.is_empty:\n            # skip when intersection is empty: this can happen for masks that consist of multiple disconnected parts\n            continue\n        if isinstance(isect, (MultiLineString, GeometryCollection)):\n            # take the line segment intersecting with `poly` that is closest to the centroid point\n            isect = isect.geoms[np.argmin([g.distance(cpoint) for g in isect.geoms])]\n        if isinstance(isect, Point):\n            # sometimes the intersection can be a singleton point\n            points.append(isect)\n            continue\n        isect = isect.boundary\n        if poly.intersects(cpoint):\n            points = [isect]\n            break\n        else:\n            points.append(isect)\n\n    if len(points) == 0:\n        # multiple reasons for this one:\n        # - if len(polys)==0\n        # - if len(polys)==1, but for some reason the line does not intersect with polygon\n        # - if the above search does not match with any points\n\n        return cx, cy\n\n    points = points[np.argmin([p.distance(cpoint) for p in points])]\n    if isinstance(points, Point):\n        return np.array(points.xy)\n    \n    points = [np.array(p.xy).tolist() for p in points.geoms]\n    visual_center = np.average(points, (0, 2))\n    return visual_center\n

\n

以下是我对这个主题的分析：

\n

\n
• get_center_of_mass是最快的，但是，正如本线程中提到的，质心可以位于非凸形状的形状之外。
\n
• get_representative_point也很快，但识别的点虽然总是保证留在形状内（或者进行较小的编辑，甚至多个断开连接的形状！），但与对象的中心没有太多关系
\n
• get_skeleton_center_of_mass返回一个感觉上不错的中心点，但速度很慢并且需要断开形状的逻辑
\n
• get_furthest_point_from_edge相对较快，很容易推广到不连续的形状，并且中心点视觉上令人愉悦
\n
• get_furthest_point_from_edge_cv其他方面的表现与以下类似，get_furthest_point_from_edge但速度快一个数量级
\n
• get_center_of_half_area_line执行得很好：结果通常最接近我自己注释视觉中心的位置。不幸的是，至少我的实现速度相当慢。
\n

\n

rows = 4\ncols = 4\nmarkers = [\'x\', \'+\', "*", "o", \'^\', "v"]\ncolors = [\'r\',\'b\',\'g\',\'orange\', \'purple\', \'lime\']\nfunctions = [\n    get_center_of_mass, \n    get_representative_point, \n    get_skeleton_center_of_mass, \n    get_furthest_point_from_edge,\n    get_furthest_point_from_edge_cv2,\n    get_center_of_half_area_line\n]\n\nplt.figure(figsize=(2*cols, 2*rows, ))\nfor i in range(rows*cols): \n    cnt = get_random_contour()\n    mask = draw_contour_on_mask((H,W), cnt)\n    \n    plt.subplot(cols,rows, i+1)\n    plt.imshow(mask, cmap=\'gray\')\n    for c, m, f in zip(colors, markers, functions):\n        l = f.__name__\n        cx, cy = f(cnt)\n        plt.scatter(cx, cy, c=c, s=100, label=l, marker=m, alpha=0.7)\n\nplt.tight_layout()    \nplt.legend(loc=3)\nplt.show()\n

\n

\n

以下是算法在 100 个随机示例上运行的速度比较：

\n

N_EXAMPLES = 100\ncnts = [get_random_contour() for _ in range(N_EXAMPLES)]\n\nfor fn in functions:\n    print(fn.__name__+":")\n    %time _ = [fn(cnt) for cnt in cnts]\n    print("~ "*40)\n

\n

get_center_of_mass:\nCPU times: user 2.35 ms, sys: 777 \xc2\xb5s, total: 3.13 ms\nWall time: 1.91 ms\n~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ \nget_representative_point:\nCPU times: user 15.7 ms, sys: 0 ns, total: 15.7 ms\nWall time: 14.8 ms\n~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ \nget_skeleton_center_of_mass:\nCPU times: user 6.52 s, sys: 104 ms, total: 6.62 s\nWall time: 6.62 s\n~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ \nget_furthest_point_from_edge:\nCPU times: user 413 ms, sys: 63 \xc2\xb5s, total: 413 ms\nWall time: 413 ms\n~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ \nget_furthest_point_from_edge_cv2:\nCPU times: user 47.8 ms, sys: 0 ns, total: 47.8 ms\nWall time: 47.8 ms\n~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ \nget_center_of_half_area_line:\nCPU times: user 1.66 s, sys: 0 ns, total: 1.66 s\nWall time: 1.66 s\n~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ \n

\n