仅考虑 3D 点，这种对 "smallest enclosing ball" 的幼稚方法是否足够好？
Considering only 3D points, wouldn't this naive approach to the "smallest enclosing ball" be good enough?
让我们在示例中使用 C#。
public class Sphere
{
public Point Center { get; set; }
public float Radius { get; set; }
public Sphere(IEnumerable<Point> points)
{
Point first = points.First();
Point vecMaxZ = first;
Point vecMinZ = first;
Point vecMaxY = first;
Point vecMinY = first;
Point vecMinX = first;
Point vecMaxX = first;
foreach (Point current in points)
{
if (current.X < vecMinX.X)
{
vecMinX = current;
}
if (current.X > vecMaxX.X)
{
vecMaxX = current;
}
if (current.Y < vecMinY.Y)
{
vecMinY = current;
}
if (current.Y > vecMaxY.Y)
{
vecMaxY = current;
}
if (current.Z < vecMinZ.Z)
{
vecMinZ = current;
}
if (current.Z > vecMaxZ.Z)
{
vecMaxZ = current;
}
}
//the lines bellow assure at least 2 points sit on the surface of the sphere.
//I'm pretty sure the algorithm is solid so far, unless I messed up the if/elses.
//I've been over this, looking at the variables and the if/elses and they all
//seem correct, but our own errors are the hardest to spot,
//so maybe there's something wrong here.
float diameterCandidateX = vecMinX.Distance(vecMaxX);
float diameterCandidateY = vecMinY.Distance(vecMaxY);
float diameterCandidateZ = vecMinZ.Distance(vecMaxZ);
Point c;
float r;
if (diameterCandidateX > diameterCandidateY)
{
if (diameterCandidateX > diameterCandidateZ)
{
c = vecMinX.Midpoint(vecMaxX);
r = diameterCandidateX / 2f;
}
else
{
c = vecMinZ.Midpoint(vecMaxZ);
r = diameterCandidateZ / 2f;
}
}
else if (diameterCandidateY > diameterCandidateZ)
{
c = vecMinY.Midpoint(vecMaxY);
r = diameterCandidateY / 2f;
}
else
{
c = vecMinZ.Midpoint(vecMaxZ);
r = diameterCandidateZ / 2f;
}
//the lines bellow look for points outside the sphere, and if one is found, then:
//1 let dist be the distance from the stray point to the current center
//2 let diff be the equal to dist - radius
//3 radius will then the increased by half of diff.
//4 a vector with the same direction as the stray point but with magnitude equal to diff is found
//5 the current center is moved by half the vector found in the step above.
//
//the stray point will now be included
//and, I would expect, the relationship between the center and other points will be mantained:
//if distance from p to center = r / k,
//then new distance from p to center' = r' / k,
//where k doesn't change from one equation to the other.
//this is where I'm wrong. I cannot figure out how to mantain this relationship.
//clearly, I'm moving the center by the wrong amount, and increasing the radius wrongly too.
//I've been over this problem for so much time, I cannot think outside the box.
//my whole world is the box. The box and I are one.
//maybe someone from outside my world (the box) could tell me where my math is wrong, please.
foreach (Point current in points)
{
float dist = current.Distance(c);
if (dist > r)
{
float diff = dist - r;
r += diff / 2f;
float scaleFactor = diff / current.Length();
Point adjust = current * scaleFactor;
c += adjust / 2f;
}
}
Center = c;
Radius = r;
}
public bool Contains(Point point) => Center.Distance(point) <= Radius;
public override string ToString() => $"Center: {Center}; Radius: {Radius}";
}
public class Point
{
public float X { get; set; }
public float Y { get; set; }
public float Z { get; set; }
public Point(float x, float y, float z)
{
X = x;
Y = y;
Z = z;
}
public float LengthSquared() => X * X + Y * Y + Z * Z;
public float Length() => (float) Math.Sqrt(X * X + Y * Y + Z * Z);
public float Distance(Point another)
{
return (float) Math.Sqrt(
(X - another.X) * (X - another.X)
+ (Y - another.Y) * (Y - another.Y)
+ (Z - another.Z) * (Z - another.Z));
}
public float DistanceSquared(Point another)
{
return (X - another.X) * (X - another.X)
+ (Y - another.Y) * (Y - another.Y)
+ (Z - another.Z) * (Z - another.Z);
}
public Point Perpendicular()
{
return new Point(-Y, X, Z);
}
public Point Midpoint(Point another)
{
return new Point(
(X + another.X) / 2f,
(Y + another.Y) / 2f,
(Z + another.Z) / 2f);
}
public override string ToString() => $"({X}, {Y}, {Z})";
public static Point operator +(Point p1, Point p2)
{
return new Point(p1.X + p2.X, p1.Y + p2.Y, p1.Z + p2.Z);
}
public static Point operator *(Point p1, float v)
{
return new Point(p1.X * v, p1.Y * v, p1.Z * v);
}
public static Point operator /(Point p1, float v)
{
return new Point(p1.X / v, p1.Y / v, p1.Z / v);
}
}
//Note: this class is here so I can be able to solve the problems suggested by
//Eric Lippert.
public class Line
{
private float coefficient;
private float constant;
public Line(Point p1, Point p2)
{
float deltaY = p2.Y - p1.Y;
float deltaX = p2.X - p1.X;
coefficient = deltaY / deltaX;
constant = coefficient * -p1.X + p1.Y;
}
public Point FromX(float x)
{
return new Point(x, x * coefficient + constant, 0);
}
public Point FromY(float y)
{
return new Point((y - constant) / coefficient, y, 0);
}
public Point Intersection(Line another)
{
float x = (another.constant - constant) / (coefficient - another.coefficient);
float y = FromX(x).Y;
return new Point(x, y, 0);
}
}
我可以安全地假设这将 运行 至少与那些通常考虑的奇特算法一样快，为了稳健性，点具有任意维数的可能性，从 2 到任何，例如 1000 或 10,000 个维度。
我只需要 3 个维度，不多也不少。由于我没有计算机科学学位（或任何与此相关的学位，我是一名高中二年级学生），因此我很难分析算法的性能和资源消耗。所以我的问题基本上是：与花哨的算法相比，我的 "smallest enclosing sphere for dumbs" 算法在性能和资源消耗方面是否良好？是否存在我的算法中断而专业算法没有中断的点，这意味着它的性能如此糟糕会导致明显的损失（例如，如果我的点数太多）。
编辑 1： 我编辑了代码，因为它根本没有意义（我很饿，当时是下午 4 点，我一整天都没有吃东西）。我认为这个更有意义，但不确定它是否正确。最初的问题是：如果这个解决了这个问题，在我们事先知道所有点都有 3 维的情况下，它是否做得足够好以与标准的专业算法竞争？
编辑 2： 现在我很确定性能很差，我对实现一个朴素的算法来找到最小的封闭球体失去了所有希望。我只想做一些有用的东西。请检查最新更新。
编辑 3: 也不起作用。我退出了。
编辑 4： 最后，我不知道……大约 5 个小时。我想到了。耶稣基督。这个有效。有人可以告诉我性能问题吗？与专业算法相比真的很糟糕吗？我可以更改哪些行以使其更好？有没有断点？请记住，我将始终将其用于 3D 点。
编辑 5： 我从 Bychenko 那里了解到，以前的算法仍然无效。我睡在这个问题上，这是我的算法的新版本。我知道它不起作用，而且我很清楚哪里错了，谁能告诉我为什么这些特定的计算是错误的以及如何解决它们？我倾向于认为这与三角函数有关。我的假设不适用于欧几里德 space，因为我无法停止将向量视为实数
一组实数，在我的例子中，我用它来确定欧几里德 space 中的一个位置。我很确定我在最后一个循环的某处遗漏了一些正弦或余弦（当然，不完全是正弦或余弦，而是笛卡尔坐标中的等价物，因为我们不知道任何角度。
EDIT 5 的附录： 关于 Eric Lippert 提出的问题：
(1) 太琐碎了：p
（2）我先为圈子做；我将为此添加 class 行。
Point a, b, c; //they are not collinear
Point midOfAB = a.Midpoint(b);
Point midOfBC = b.Midpoint(c);
//multiplying the vector by a scalar as I do bellow doesn't matter right?
Point perpendicularToAB = midOfAB.Perpendicular() * 3;
Point perpendicularToBC = midOfBC.Perpendicular() * 3;
Line bisectorAB = new Line(perpendicularToAB, midOfAB);
Line bisectorBC = new Line(perpendicularToBC, midOfBC);
Point center = bisectorAB.Intersection(bisectorBC);
float distA = center.Distance(a);
float distB = center.Distance(b);
float distC = center.Distance(c);
if(distA == distB && distB == distC)
//it works (spoiler alert: it doesn't)
else
//you're a failure, programmer, pick up your skate and practice some ollies
抱歉，您的算法错误。没有解决问题。
反例（3分）：
A = (0, 0, 0) - closest to origin (0)
B = (3, 3, 0) - farthest from origin (3 * sqrt(2) == 4.2426...)
C = (4, 0, 0)
您的天真算法声明球体的中心位于
P = (3 / sqrt(2), 3 / sqrt(2), 0)
和半径
R = 3 / sqrt(2)
你可以看到点 C = (4, 0, 0) 超出 球体
编辑更新的（但天真的）算法仍然是错误的。
反例（3分）：
A = (0, 0, 0)
B = (1, 2, 0)
C = (4, 1, 0)
根据算法，球体的中心位于
P = (2, 1, 0)
半径
R = sqrt(5)
你可以看到球体不是 minimal（最小的）。
第 N 次编辑 你的算法仍然不正确。在探索 灰色地带 （您知道问题所在，但部分存在漏洞）时，投资测试自动化是一个很好的做法。如您所知，在三角形的情况下，所有顶点都应在球体上；让我们根据这个事实验证您的解决方案：
public static class SphereValidator {
private static Random m_Random = new Random();
private static String Validate() {
var triangle = Enumerable
.Range(0, 3)
.Select(i => new Point(m_Random.Next(100), m_Random.Next(100), m_Random.Next(100)))
.ToArray();
Sphere solution = new Sphere(triangle);
double tolerance = 1.0e-5;
for (int i = 0; i < triangle.Length; ++i) {
double r = triangle[i].Distance(solution.Center);
if (Math.Abs(r - solution.Radius) > tolerance) {
return String.Format("Counter example\r\n A: {0}\r\n B: {1}\r\n C: {2}\r\n expected distance to \"{3}\": {4}; actual R {5}",
triangle[0], triangle[1], triangle[2], (char) ('A' + i), r, solution.Radius);
}
}
return null;
}
public static String FindCounterExample(int attempts = 10000) {
for (int i = 0; i < attempts; ++i) {
String result = Validate();
if (!String.IsNullOrEmpty(result))
Console.WriteLine(result);
return;
}
Console.WriteLine(String.Format("Yes! All {0} tests passed!", attempts));
}
}
我刚刚 运行 上面的代码并得到：
Counter example
A: (3, 30, 9)
B: (1, 63, 40)
C: (69, 1, 16)
expected distance to "A": 35.120849609375; actual R 53.62698
对于粗略的近似，计算轴对齐边界框，然后计算该框的边界球（同心，直径 = √(W² + H² + D²)）。
您可以通过计算与该中心的最大距离来进行优化。
让我们在示例中使用 C#。
public class Sphere
{
public Point Center { get; set; }
public float Radius { get; set; }
public Sphere(IEnumerable<Point> points)
{
Point first = points.First();
Point vecMaxZ = first;
Point vecMinZ = first;
Point vecMaxY = first;
Point vecMinY = first;
Point vecMinX = first;
Point vecMaxX = first;
foreach (Point current in points)
{
if (current.X < vecMinX.X)
{
vecMinX = current;
}
if (current.X > vecMaxX.X)
{
vecMaxX = current;
}
if (current.Y < vecMinY.Y)
{
vecMinY = current;
}
if (current.Y > vecMaxY.Y)
{
vecMaxY = current;
}
if (current.Z < vecMinZ.Z)
{
vecMinZ = current;
}
if (current.Z > vecMaxZ.Z)
{
vecMaxZ = current;
}
}
//the lines bellow assure at least 2 points sit on the surface of the sphere.
//I'm pretty sure the algorithm is solid so far, unless I messed up the if/elses.
//I've been over this, looking at the variables and the if/elses and they all
//seem correct, but our own errors are the hardest to spot,
//so maybe there's something wrong here.
float diameterCandidateX = vecMinX.Distance(vecMaxX);
float diameterCandidateY = vecMinY.Distance(vecMaxY);
float diameterCandidateZ = vecMinZ.Distance(vecMaxZ);
Point c;
float r;
if (diameterCandidateX > diameterCandidateY)
{
if (diameterCandidateX > diameterCandidateZ)
{
c = vecMinX.Midpoint(vecMaxX);
r = diameterCandidateX / 2f;
}
else
{
c = vecMinZ.Midpoint(vecMaxZ);
r = diameterCandidateZ / 2f;
}
}
else if (diameterCandidateY > diameterCandidateZ)
{
c = vecMinY.Midpoint(vecMaxY);
r = diameterCandidateY / 2f;
}
else
{
c = vecMinZ.Midpoint(vecMaxZ);
r = diameterCandidateZ / 2f;
}
//the lines bellow look for points outside the sphere, and if one is found, then:
//1 let dist be the distance from the stray point to the current center
//2 let diff be the equal to dist - radius
//3 radius will then the increased by half of diff.
//4 a vector with the same direction as the stray point but with magnitude equal to diff is found
//5 the current center is moved by half the vector found in the step above.
//
//the stray point will now be included
//and, I would expect, the relationship between the center and other points will be mantained:
//if distance from p to center = r / k,
//then new distance from p to center' = r' / k,
//where k doesn't change from one equation to the other.
//this is where I'm wrong. I cannot figure out how to mantain this relationship.
//clearly, I'm moving the center by the wrong amount, and increasing the radius wrongly too.
//I've been over this problem for so much time, I cannot think outside the box.
//my whole world is the box. The box and I are one.
//maybe someone from outside my world (the box) could tell me where my math is wrong, please.
foreach (Point current in points)
{
float dist = current.Distance(c);
if (dist > r)
{
float diff = dist - r;
r += diff / 2f;
float scaleFactor = diff / current.Length();
Point adjust = current * scaleFactor;
c += adjust / 2f;
}
}
Center = c;
Radius = r;
}
public bool Contains(Point point) => Center.Distance(point) <= Radius;
public override string ToString() => $"Center: {Center}; Radius: {Radius}";
}
public class Point
{
public float X { get; set; }
public float Y { get; set; }
public float Z { get; set; }
public Point(float x, float y, float z)
{
X = x;
Y = y;
Z = z;
}
public float LengthSquared() => X * X + Y * Y + Z * Z;
public float Length() => (float) Math.Sqrt(X * X + Y * Y + Z * Z);
public float Distance(Point another)
{
return (float) Math.Sqrt(
(X - another.X) * (X - another.X)
+ (Y - another.Y) * (Y - another.Y)
+ (Z - another.Z) * (Z - another.Z));
}
public float DistanceSquared(Point another)
{
return (X - another.X) * (X - another.X)
+ (Y - another.Y) * (Y - another.Y)
+ (Z - another.Z) * (Z - another.Z);
}
public Point Perpendicular()
{
return new Point(-Y, X, Z);
}
public Point Midpoint(Point another)
{
return new Point(
(X + another.X) / 2f,
(Y + another.Y) / 2f,
(Z + another.Z) / 2f);
}
public override string ToString() => $"({X}, {Y}, {Z})";
public static Point operator +(Point p1, Point p2)
{
return new Point(p1.X + p2.X, p1.Y + p2.Y, p1.Z + p2.Z);
}
public static Point operator *(Point p1, float v)
{
return new Point(p1.X * v, p1.Y * v, p1.Z * v);
}
public static Point operator /(Point p1, float v)
{
return new Point(p1.X / v, p1.Y / v, p1.Z / v);
}
}
//Note: this class is here so I can be able to solve the problems suggested by
//Eric Lippert.
public class Line
{
private float coefficient;
private float constant;
public Line(Point p1, Point p2)
{
float deltaY = p2.Y - p1.Y;
float deltaX = p2.X - p1.X;
coefficient = deltaY / deltaX;
constant = coefficient * -p1.X + p1.Y;
}
public Point FromX(float x)
{
return new Point(x, x * coefficient + constant, 0);
}
public Point FromY(float y)
{
return new Point((y - constant) / coefficient, y, 0);
}
public Point Intersection(Line another)
{
float x = (another.constant - constant) / (coefficient - another.coefficient);
float y = FromX(x).Y;
return new Point(x, y, 0);
}
}
我可以安全地假设这将 运行 至少与那些通常考虑的奇特算法一样快，为了稳健性，点具有任意维数的可能性，从 2 到任何，例如 1000 或 10,000 个维度。
我只需要 3 个维度，不多也不少。由于我没有计算机科学学位（或任何与此相关的学位，我是一名高中二年级学生），因此我很难分析算法的性能和资源消耗。所以我的问题基本上是：与花哨的算法相比，我的 "smallest enclosing sphere for dumbs" 算法在性能和资源消耗方面是否良好？是否存在我的算法中断而专业算法没有中断的点，这意味着它的性能如此糟糕会导致明显的损失（例如，如果我的点数太多）。
编辑 1： 我编辑了代码，因为它根本没有意义（我很饿，当时是下午 4 点，我一整天都没有吃东西）。我认为这个更有意义，但不确定它是否正确。最初的问题是：如果这个解决了这个问题，在我们事先知道所有点都有 3 维的情况下，它是否做得足够好以与标准的专业算法竞争？
编辑 2： 现在我很确定性能很差，我对实现一个朴素的算法来找到最小的封闭球体失去了所有希望。我只想做一些有用的东西。请检查最新更新。
编辑 3: 也不起作用。我退出了。
编辑 4： 最后，我不知道……大约 5 个小时。我想到了。耶稣基督。这个有效。有人可以告诉我性能问题吗？与专业算法相比真的很糟糕吗？我可以更改哪些行以使其更好？有没有断点？请记住，我将始终将其用于 3D 点。
编辑 5： 我从 Bychenko 那里了解到，以前的算法仍然无效。我睡在这个问题上，这是我的算法的新版本。我知道它不起作用，而且我很清楚哪里错了，谁能告诉我为什么这些特定的计算是错误的以及如何解决它们？我倾向于认为这与三角函数有关。我的假设不适用于欧几里德 space，因为我无法停止将向量视为实数 一组实数，在我的例子中，我用它来确定欧几里德 space 中的一个位置。我很确定我在最后一个循环的某处遗漏了一些正弦或余弦（当然，不完全是正弦或余弦，而是笛卡尔坐标中的等价物，因为我们不知道任何角度。
EDIT 5 的附录： 关于 Eric Lippert 提出的问题： (1) 太琐碎了：p （2）我先为圈子做；我将为此添加 class 行。
Point a, b, c; //they are not collinear
Point midOfAB = a.Midpoint(b);
Point midOfBC = b.Midpoint(c);
//multiplying the vector by a scalar as I do bellow doesn't matter right?
Point perpendicularToAB = midOfAB.Perpendicular() * 3;
Point perpendicularToBC = midOfBC.Perpendicular() * 3;
Line bisectorAB = new Line(perpendicularToAB, midOfAB);
Line bisectorBC = new Line(perpendicularToBC, midOfBC);
Point center = bisectorAB.Intersection(bisectorBC);
float distA = center.Distance(a);
float distB = center.Distance(b);
float distC = center.Distance(c);
if(distA == distB && distB == distC)
//it works (spoiler alert: it doesn't)
else
//you're a failure, programmer, pick up your skate and practice some ollies
抱歉，您的算法错误。没有解决问题。
反例（3分）：
A = (0, 0, 0) - closest to origin (0)
B = (3, 3, 0) - farthest from origin (3 * sqrt(2) == 4.2426...)
C = (4, 0, 0)
您的天真算法声明球体的中心位于
P = (3 / sqrt(2), 3 / sqrt(2), 0)
和半径
R = 3 / sqrt(2)
你可以看到点 C = (4, 0, 0) 超出 球体
编辑更新的（但天真的）算法仍然是错误的。
反例（3分）：
A = (0, 0, 0)
B = (1, 2, 0)
C = (4, 1, 0)
根据算法，球体的中心位于
P = (2, 1, 0)
半径
R = sqrt(5)
你可以看到球体不是 minimal（最小的）。
第 N 次编辑 你的算法仍然不正确。在探索 灰色地带 （您知道问题所在，但部分存在漏洞）时，投资测试自动化是一个很好的做法。如您所知，在三角形的情况下，所有顶点都应在球体上；让我们根据这个事实验证您的解决方案：
public static class SphereValidator {
private static Random m_Random = new Random();
private static String Validate() {
var triangle = Enumerable
.Range(0, 3)
.Select(i => new Point(m_Random.Next(100), m_Random.Next(100), m_Random.Next(100)))
.ToArray();
Sphere solution = new Sphere(triangle);
double tolerance = 1.0e-5;
for (int i = 0; i < triangle.Length; ++i) {
double r = triangle[i].Distance(solution.Center);
if (Math.Abs(r - solution.Radius) > tolerance) {
return String.Format("Counter example\r\n A: {0}\r\n B: {1}\r\n C: {2}\r\n expected distance to \"{3}\": {4}; actual R {5}",
triangle[0], triangle[1], triangle[2], (char) ('A' + i), r, solution.Radius);
}
}
return null;
}
public static String FindCounterExample(int attempts = 10000) {
for (int i = 0; i < attempts; ++i) {
String result = Validate();
if (!String.IsNullOrEmpty(result))
Console.WriteLine(result);
return;
}
Console.WriteLine(String.Format("Yes! All {0} tests passed!", attempts));
}
}
我刚刚 运行 上面的代码并得到：
Counter example
A: (3, 30, 9)
B: (1, 63, 40)
C: (69, 1, 16)
expected distance to "A": 35.120849609375; actual R 53.62698
对于粗略的近似，计算轴对齐边界框，然后计算该框的边界球（同心，直径 = √(W² + H² + D²)）。
您可以通过计算与该中心的最大距离来进行优化。