已知空间上若干点(xi, yi, zi), 求空间上包含这些点的最小球半径 R, 以及球心坐标。
思路:球心与这些点的最大距离为半径, 球心与最大距离点生成向量,将球心朝着该向量方向移动若干距离,再计算半径的变化。
namespace Test_BST
{
public class Program
{
static void Main(string[] args)
{
// 初始化输入点
List<Point> originPoints = new List<Point>() { ............};
double radius = AnnealAlgorithm(originPoints);
}
private struct Point
{
public double x;
public double y;
public double z;
}
// square of a number
private static double Sqr(double x) { return x * x; }
// 两点之间的距离
private static double Dist(Point A, Point B)
{
return Math.Sqrt(Sqr(A.x - B.x) + Sqr(A.y - B.y) + Sqr(A.z - B.z));
}
// 求最大半径
private static double GetMaxRadius(Point p0, List<Point> pts)
{
double maxRadius = 0;
foreach (var point in pts)
{
double radius = Dist(p0, point);
maxRadius = radius > maxRadius ? radius : maxRadius;
}
return maxRadius;
}
private static double AnnealAlgorithm(List<Point> originPts)
{
Point center = new Point();
center.x = 0;
center.y = 0;
center.z = 0;
// 将初始化中心点设置为所有点的代数平均位置
foreach (var pt in originPts)
{
center.x += pt.x;
center.y += pt.y;
center.z += pt.z;
}
center.x /= originPts.Count;
center.y /= originPts.Count;
center.z /= originPts.Count;
double temp = 1e3; // 初始温度
double coolingFactor = 0.98; // 降温因子
double ans = GetMaxRadius(center, originPts); // 当前最小半径
var random = new Random();
while (temp > 1e-5)
{
Point newCenter = new Point();
double max_r = 0;
// 找到与当前中心点距离最远的点,将中心向着改点移动
for (int i = 0; i < originPts.Count; i++)
{
double r = Dist(center, originPts[i]);
if (r > max_r)
{
newCenter.x = (originPts[i].x - center.x) / r;
newCenter.y = (originPts[i].y - center.y) / r;
newCenter.z = (originPts[i].z - center.z) / r;
max_r = r;
}
}
newCenter.x = center.x + newCenter.x * temp;
newCenter.y = center.y + newCenter.y * temp;
newCenter.z = center.z + newCenter.z * temp;
// 移动后的最大半径
double tmp = GetMaxRadius(newCenter, originPts);
if (tmp < ans)
{
center.x += newCenter.x * temp;
center.y += newCenter.y * temp;
center.z += newCenter.z * temp;
}
else if (Math.Exp((ans -tmp)/temp) > random.NextDouble() )
{
center.x += newCenter.x * temp;
center.y += newCenter.y * temp;
center.z += newCenter.z * temp;
}
temp *= coolingFactor;
}
double miniRadius = GetMaxRadius(center, originPts);
Console.WriteLine("the cooridnate of the center is {0}, the radius value is {1}", center, miniRadius));
return miniRadius;
}
}
}