Python中最小的圆圈，代码错误 [英] Smallest enclosing circle in Python, error in the code
问题描述
我有一组表示多边形的顶点（x，y）的点。
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)]
我需要找到最小的封闭圆（中心的面积，x和y，和半径）
I need to find the Smallest Enclosing Circle (area, x and y of center, and radius)
我正在使用从此页面派生的python代码： Nayuki的最小包围圈
I am using the python code derived from this page: Smallest enclosing circle of Nayuki
当我运行代码时，结果会在每个ti我，例如：
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)
其中x是y （圆心）和半径
where return is x, y (of the center of the circle), and radius
使用商业软件（例如ArcGiS）在一些点集上得出正确的结果是：
using a commercial software (i.e. ArcGiS) whin the some set of points the correct result is:
421646.74552, 4596389.82475, 24.323246
使用的代码：
#
# 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)
推荐答案
在不了解您的算法的情况下，我注意到了一件事：坐标之间的比率并且您的半径很大，大约2e5。也许，当试图在离原点如此远的点周围找到一个圆时，您的算法不适。尤其是在您的 _make_circumcircle 函数中，这会导致大数减法，这对于数字错误通常是一件坏事。
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])
结果：
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).
这篇关于Python中最小的圆圈，代码错误的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
