最小的封闭圆，代码错误 [英] Smallest enclosing circle, error in the code

问题描述

我有一组表示多边形顶点 (x, y) 的点.

<预> <代码>点= [(421640.3639270504，4596366.353552659)，(421635.79361391126，4596369.054192241)，(421632.6774913164，4596371.131607305)，(421629.14588570886，4596374.870954419)，(421625.6142801013，4596377.779335507)，(421624.99105558236，4596382.14190714)，(421630.1845932406，4596388.062540068)，(421633.3007158355，4596388.270281575)，(421637.87102897465，4596391.8018871825)，(421642.4413421138，4596394.918009778)，(421646.5961722403，4596399.903805929)，(421649.71229483513，4596403.850894549)，(421653.8940752105，4596409.600842565)，(421654.69809098693，4596410.706364258)，(421657.60647207545，4596411.329588776)，(421660.514853164，4596409.875398233)，(421661.3458191893，4596406.136051118)，(421661.5535606956，4596403.22767003)，(421658.85292111343，4596400.94251346)，(421656.5677645438，4596399.696064423)，(421655.52905701223，4596396.164458815)，(421652.82841743，4596394.502526765)，(421648.46584579715，4596391.8018871825)，(421646.38843073393，4596388.270281575)、(421645.55746)470863，4596386.400608018)，(421647.21939675923，4596384.115451449)，(421649.5045533288，4596382.661260904)，(421650.7510023668，4596378.714172284)，(421647.8426212782，4596375.8057911955)，(421644.9342401897，4596372.897410107)，(421643.6877911517，4596370.404512031)，(421640.3639270504，4596366.353552659)]

我需要找到最小的封闭圆(面积、中心的 x 和 y 以及半径)

我正在使用从这个页面派生的 python 代码:最小的封闭圈奈雪

当我运行代码时，结果每次都会改变，例如:

<预><代码>>>>make_circle(点)(421643.0645666326, 4596393.82736687, 23.70763190712525)>>>make_circle(点)(421647.8426212782, 4596375.8057911955, 0.0)>>>make_circle(点)(421648.9851995629, 4596388.841570718, 24.083963460031157)

其中 return 是 x、y(圆心)和半径

使用商业软件(即 ArcGiS)时，正确的结果是:

421646.74552、4596389.82475、24.323246

使用的代码:

<预><代码>## 最小的封闭圆## 版权所有 (c) 2014 Project Nayuki# https://www.nayuki.io/page/smallest-enclosure-circle## 该程序是免费软件:您可以重新分发和/或修改它# 它根据由发布的 GNU 通用公共许可证的条款# 自由软件基金会，许可证的第 3 版，或#(由您选择)任何更高版本.## 本程序已分发，希望对大家有用，# 但没有任何保证；甚至没有暗示的保证# 特定用途的适销性或适用性.见# GNU 通用公共许可证了解更多详情.## 您应该已经收到一份 GNU 通用公共许可证的副本# 连同这个程序(见COPYING.txt).# 如果没有，请参阅 <http://www.gnu.org/licenses/>.#导入数学，随机# 数据约定:一个点是一对浮点数 (x, y).圆是三个浮点数(中心 x，中心 y，半径).## 返回包含所有给定点的最小圆.在预期的 O(n) 时间内运行，随机.# 输入:一系列浮点数或整数对，例如[(0,5), (3.1,-2.7)].# 输出:代表一个圆的三个浮点数.# 注意:如果给0分，则返回None.如果给定 1 点，则返回半径为 0 的圆.#def make_circle(点):# 转换为浮动和随机顺序shuffled = [(float(p[0]), float(p[1])) for p in points]random.shuffle(shuffled)# 逐步向圆添加点或重新计算圆c = 无for (i, p) in enumerate(shuffled):如果 c 是 None 或不是 _is_in_circle(c, p):c = _make_circle_one_point(shuffled[0: i + 1], p)返回 c# 已知一个边界点def _make_circle_one_point(points, p):c = (p[0], p[1], 0.0)对于枚举(点)中的(i，q):如果不是 _is_in_circle(c, q):如果 c[2] == 0.0:c = _make_diameter(p, q)别的:c = _make_circle_two_points(points[0: i + 1], p, q)返回 c# 已知两个边界点def _make_circle_two_points(points, p, q):直径 = _make_diameter(p, q)if all(_is_in_circle(diameter, r) for r in points):回程直径左 = 无正确 = 无对于 r 点:cross = _cross_product(p[0], p[1], q[0], q[1], r[0], r[1])c = _make_circumcircle(p, q, r)如果 c 是 None:继续elif 交叉 >0.0 and (left is None or _cross_product(p[0], p[1], q[0], q[1], c[0], c[1]) > _cross_product(p[0], p[1], q[0], q[1], left[0], left[1])):左 = celif 交叉 <0.0 and (right is None or _cross_product(p[0], p[1], q[0], q[1], c[0], c[1]) <_cross_product(p[0], p[1], q[0], q[1], right[0], right[1])):右 = creturn left if (right is None or (left is not None and left[2] <= right[2])) else rightdef_make_circumcircle(p0, p1, p2):# 维基百科的数学算法:外接圆ax = p0[0];ay = p0[1]bx = p1[0];通过 = p1[1]cx = p2[0];cy = p2[1]d = (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)) * 2.0如果 d == 0.0:返回无x = ((ax * ax + ay * ay) * (by - cy) + (bx * bx + by * by) * (cy - ay) + (cx * cx + cy * cy) * (ay - by))/dy = ((ax * ax + ay * ay) * (cx - bx) + (bx * bx + by * by) * (ax - cx) + (cx * cx + cy * cy) * (bx - ax))/d返回 (x, y, math.hypot(x - ax, y - ay))def _make_diameter(p0, p1):返回 ((p0[0] + p1[0])/2.0, (p0[1] + p1[1])/2.0, math.hypot(p0[0] - p1[0], p0[1] - p1[1])/2.0)_EPSILON = 1e-12def _is_in_circle(c, p):返回 c 不是 None 并且 math.hypot(p[0] - c[0], p[1] - c[1])

解决方案

在不了解您的算法的情况下，我注意到一件事:您的坐标和半径之间的比率非常大，大约为 2e5.也许，您的算法在尝试在远离原点的点周围找到一个圆时条件不佳.尤其是在您的 _make_circumcircle 函数中，这会导致大数相减，这对于数值错误通常是一件坏事.

由于相对于点拟合半径和圆心应该与平移无关，您可以简单地减去所有点的平均值(点云的质心)，进行拟合，然后将平均值相加得到最终结果:

def numeric_stable_circle(points):pts = np.array(点)mean_pts = np.mean(pts, 0)打印'点的平均值:'，mean_ptspts -= mean_pts # 向原点平移结果 = make_circle(pts)打印 'result without mean:', result打印 'result with mean:', (result[0] + mean_pts[0],结果[1] + mean_pts[1], 结果[2])

结果:

平均分:[421645.83745955 4596388.99204294]结果无均值:(0.9080813432488977, 0.8327111343034483, 24.323287017466253)结果均值:(421646.74554089626, 4596389.8247540779, 24.323287017466253)

这些数字不会从一次运行到下一次改变一个数字，并且与您的正确结果"只有很小的差异(可能由于不同的实现而导致不同的数字错误).

I have a set of points that represent the vertices (x, y) of a polygon.

points= [(421640.3639270504, 4596366.353552659), (421635.79361391126, 4596369.054192241), (421632.6774913164, 4596371.131607305), (421629.14588570886, 4596374.870954419), (421625.6142801013, 4596377.779335507), (421624.99105558236, 4596382.14190714), (421630.1845932406, 4596388.062540068), (421633.3007158355, 4596388.270281575), (421637.87102897465, 4596391.8018871825), (421642.4413421138, 4596394.918009778), (421646.5961722403, 4596399.903805929), (421649.71229483513, 4596403.850894549), (421653.8940752105, 4596409.600842565), (421654.69809098693, 4596410.706364258), (421657.60647207545, 4596411.329588776), (421660.514853164, 4596409.875398233), (421661.3458191893, 4596406.136051118), (421661.5535606956, 4596403.22767003), (421658.85292111343, 4596400.94251346), (421656.5677645438, 4596399.696064423), (421655.52905701223, 4596396.164458815), (421652.82841743, 4596394.502526765), (421648.46584579715, 4596391.8018871825), (421646.38843073393, 4596388.270281575), (421645.55746470863, 4596386.400608018), (421647.21939675923, 4596384.115451449), (421649.5045533288, 4596382.661260904), (421650.7510023668, 4596378.714172284), (421647.8426212782, 4596375.8057911955), (421644.9342401897, 4596372.897410107), (421643.6877911517, 4596370.404512031), (421640.3639270504, 4596366.353552659)]

I need to find the Smallest Enclosing Circle (area, x and y of center, and radius)

I am using the python code derived from this page: Smallest enclosing circle of Nayuki

when I run the code the results change every time, for example:

>>> make_circle(points)
(421643.0645666326, 4596393.82736687, 23.70763190712525)
>>> make_circle(points)
(421647.8426212782, 4596375.8057911955, 0.0)
>>> make_circle(points)
(421648.9851995629, 4596388.841570718, 24.083963460031157)

where return is x, y (of the center of the circle), and radius

using a commercial software (i.e. ArcGiS) whin the some set of points the correct result is:

421646.74552, 4596389.82475, 24.323246

code used:

# 
# Smallest enclosing circle
# 
# Copyright (c) 2014 Project Nayuki
# https://www.nayuki.io/page/smallest-enclosing-circle
# 
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program (see COPYING.txt).
# If not, see <http://www.gnu.org/licenses/>.
# 

import math, random


# Data conventions: A point is a pair of floats (x, y). A circle is a triple of floats (center x, center y, radius).

# 
# Returns the smallest circle that encloses all the given points. Runs in expected O(n) time, randomized.
# Input: A sequence of pairs of floats or ints, e.g. [(0,5), (3.1,-2.7)].
# Output: A triple of floats representing a circle.
# Note: If 0 points are given, None is returned. If 1 point is given, a circle of radius 0 is returned.
# 
def make_circle(points):
    # Convert to float and randomize order
    shuffled = [(float(p[0]), float(p[1])) for p in points]
    random.shuffle(shuffled)

    # Progressively add points to circle or recompute circle
    c = None
    for (i, p) in enumerate(shuffled):
        if c is None or not _is_in_circle(c, p):
            c = _make_circle_one_point(shuffled[0 : i + 1], p)
    return c


# One boundary point known
def _make_circle_one_point(points, p):
    c = (p[0], p[1], 0.0)
    for (i, q) in enumerate(points):
        if not _is_in_circle(c, q):
            if c[2] == 0.0:
                c = _make_diameter(p, q)
            else:
                c = _make_circle_two_points(points[0 : i + 1], p, q)
    return c


# Two boundary points known
def _make_circle_two_points(points, p, q):
    diameter = _make_diameter(p, q)
    if all(_is_in_circle(diameter, r) for r in points):
        return diameter

    left = None
    right = None
    for r in points:
        cross = _cross_product(p[0], p[1], q[0], q[1], r[0], r[1])
        c = _make_circumcircle(p, q, r)
        if c is None:
            continue
        elif cross > 0.0 and (left is None or _cross_product(p[0], p[1], q[0], q[1], c[0], c[1]) > _cross_product(p[0], p[1], q[0], q[1], left[0], left[1])):
            left = c
        elif cross < 0.0 and (right is None or _cross_product(p[0], p[1], q[0], q[1], c[0], c[1]) < _cross_product(p[0], p[1], q[0], q[1], right[0], right[1])):
            right = c
    return left if (right is None or (left is not None and left[2] <= right[2])) else right


def _make_circumcircle(p0, p1, p2):
    # Mathematical algorithm from Wikipedia: Circumscribed circle
    ax = p0[0]; ay = p0[1]
    bx = p1[0]; by = p1[1]
    cx = p2[0]; cy = p2[1]
    d = (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)) * 2.0
    if d == 0.0:
        return None
    x = ((ax * ax + ay * ay) * (by - cy) + (bx * bx + by * by) * (cy - ay) + (cx * cx + cy * cy) * (ay - by)) / d
    y = ((ax * ax + ay * ay) * (cx - bx) + (bx * bx + by * by) * (ax - cx) + (cx * cx + cy * cy) * (bx - ax)) / d
    return (x, y, math.hypot(x - ax, y - ay))


def _make_diameter(p0, p1):
    return ((p0[0] + p1[0]) / 2.0, (p0[1] + p1[1]) / 2.0, math.hypot(p0[0] - p1[0], p0[1] - p1[1]) / 2.0)


_EPSILON = 1e-12

def _is_in_circle(c, p):
    return c is not None and math.hypot(p[0] - c[0], p[1] - c[1]) < c[2] + _EPSILON


# Returns twice the signed area of the triangle defined by (x0, y0), (x1, y1), (x2, y2)
def _cross_product(x0, y0, x1, y1, x2, y2):
    return (x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0)

解决方案

Without understanding anything about your algorithm, I noticed one thing: the ratio between your coordinates and your radius is very large, about 2e5. Maybe, your algorithm is ill conditioned when trying to find a circle around points which are so far away from the origin. Especially in your _make_circumcircle function, this leads to the subtraction of large numbers, which is usually a bad thing for numerical errors.

Since fitting the radius and the center of the circle with respect to the points should be independent of a translation, you could simply subtract the mean of all points (the center of mass of your cloud of points), do the fitting, and then add the mean back to obtain the final result:

def numerical_stable_circle(points):
    pts = np.array(points)
    mean_pts = np.mean(pts, 0)
    print 'mean of points:', mean_pts
    pts -= mean_pts  # translate towards origin
    result = make_circle(pts)
    print 'result without mean:', result
    print 'result with mean:', (result[0] + mean_pts[0], 
    result[1] + mean_pts[1], result[2])

Result:

mean of points: [  421645.83745955  4596388.99204294]
result without mean: (0.9080813432488977, 0.8327111343034483, 24.323287017466253)
result with mean: (421646.74554089626, 4596389.8247540779, 24.323287017466253)

These numbers do not change a single digit from one run to the next one, and differ from your 'correct result' by only a tiny amount (probably different numerical errors due to a different implementation).

这篇关于最小的封闭圆，代码错误的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！