问题是我对编程非常陌生，现在我需要制作一个程序来计算许多球体的中心（最多 36 个，最小 3 个），每个球体有 4 个点 X，Y，Z。因为我的程序读取了一个 TXT 文件，其中包含我将其存储在列表中的点数据，其结构如下
bolas[n].xyz[row,element]
这意味着我的球体 1 的第一组点如下所示：
bolas[0] =
row0. -> [0] [1] [2]
row1. -> [0] [1] [2]
row2. -> [0] [1] [2]
row3. -> [0] [1] [2]
因此，如果我想使用球体中第 1 行中的 X 值，我必须这样做：
bolas[0].xyz[0,0]
在网上搜索时，我发现有人转换java代码并将其实现为c#来计算球体的中心，他创建了一个类，但我很新，我不知道如何使用他的类上的元素，我应该如何将我的数据引入他的班级以及如何获得结果；这是课程：
/// <summary>
/// Given four points in 3D space, solves for a sphere such that all four points
/// lie on the sphere's surface.
/// </summary>
/// <remarks>
/// Translated from Javascript on http://www.convertalot.com/sphere_solver.html, originally
/// linked to by http://stackoverflow.com/questions/13600739/calculate-centre-of-sphere-whose-surface-contains-4-points-c.
/// </remarks>
public class CircumcentreSolver
{
private const float ZERO = 0;
private double m_X0, m_Y0, m_Z0;
private double m_Radius;
private double[,] P =
{
{ ZERO, ZERO, ZERO },
{ ZERO, ZERO, ZERO },
{ ZERO, ZERO, ZERO },
{ ZERO, ZERO, ZERO }
};
/// <summary>
/// The centre of the resulting sphere.
/// </summary>
public double[] Centre
{
get { return new double[] { this.m_X0, this.m_Y0, this.m_Z0 }; }
}
/// <summary>
/// The radius of the resulting sphere.
/// </summary>
public double Radius
{
get { return this.m_Radius; }
}
/// <summary>
/// Whether the result was a valid sphere.
/// </summary>
public bool Valid
{
get { return this.m_Radius != 0; }
}
/// <summary>
/// Computes the centre of a sphere such that all four specified points in
/// 3D space lie on the sphere's surface.
/// </summary>
/// <param name="a">The first point (array of 3 doubles for X, Y, Z).</param>
/// <param name="b">The second point (array of 3 doubles for X, Y, Z).</param>
/// <param name="c">The third point (array of 3 doubles for X, Y, Z).</param>
/// <param name="d">The fourth point (array of 3 doubles for X, Y, Z).</param>
public CircumcentreSolver(double[] a, double[] b, double[] c, double[] d)
{
this.Compute(a, b, c, d);
}
/// <summary>
/// Evaluate the determinant.
/// </summary>
private void Compute(double[] a, double[] b, double[] c, double[] d)
{
P[0, 0] = a[0];
P[0, 1] = a[1];
P[0, 2] = a[2];
P[1, 0] = b[0];
P[1, 1] = b[1];
P[1, 2] = b[2];
P[2, 0] = c[0];
P[2, 1] = c[1];
P[2, 2] = c[2];
P[3, 0] = d[0];
P[3, 1] = d[1];
P[3, 2] = d[2];
// Compute result sphere.
this.Sphere();
}
private void Sphere()
{
double r, m11, m12, m13, m14, m15;
double[,] a =
{
{ ZERO, ZERO, ZERO, ZERO },
{ ZERO, ZERO, ZERO, ZERO },
{ ZERO, ZERO, ZERO, ZERO },
{ ZERO, ZERO, ZERO, ZERO }
};
// Find minor 1, 1.
for (int i = 0; i < 4; i++)
{
a[i, 0] = P[i, 0];
a[i, 1] = P[i, 1];
a[i, 2] = P[i, 2];
a[i, 3] = 1;
}
m11 = this.Determinant(a, 4);
// Find minor 1, 2.
for (int i = 0; i < 4; i++)
{
a[i, 0] = P[i, 0] * P[i, 0] + P[i, 1] * P[i, 1] + P[i, 2] * P[i, 2];
a[i, 1] = P[i, 1];
a[i, 2] = P[i, 2];
a[i, 3] = 1;
}
m12 = this.Determinant(a, 4);
// Find minor 1, 3.
for (int i = 0; i < 4; i++)
{
a[i, 0] = P[i, 0] * P[i, 0] + P[i, 1] * P[i, 1] + P[i, 2] * P[i, 2];
a[i, 1] = P[i, 0];
a[i, 2] = P[i, 2];
a[i, 3] = 1;
}
m13 = this.Determinant(a, 4);
// Find minor 1, 4.
for (int i = 0; i < 4; i++)
{
a[i, 0] = P[i, 0] * P[i, 0] + P[i, 1] * P[i, 1] + P[i, 2] * P[i, 2];
a[i, 1] = P[i, 0];
a[i, 2] = P[i, 1];
a[i, 3] = 1;
}
m14 = this.Determinant(a, 4);
// Find minor 1, 5.
for (int i = 0; i < 4; i++)
{
a[i, 0] = P[i, 0] * P[i, 0] + P[i, 1] * P[i, 1] + P[i, 2] * P[i, 2];
a[i, 1] = P[i, 0];
a[i, 2] = P[i, 1];
a[i, 3] = P[i, 2];
}
m15 = this.Determinant(a, 4);
// Calculate result.
if (m11 == 0)
{
this.m_X0 = 0;
this.m_Y0 = 0;
this.m_Z0 = 0;
this.m_Radius = 0;
}
else
{
this.m_X0 = 0.5 * m12 / m11;
this.m_Y0 = -0.5 * m13 / m11;
this.m_Z0 = 0.5 * m14 / m11;
this.m_Radius = System.Math.Sqrt(this.m_X0 * this.m_X0 + this.m_Y0 * this.m_Y0 + this.m_Z0 * this.m_Z0 - m15 / m11);
}
}
/// <summary>
/// Recursive definition of determinate using expansion by minors.
/// </summary>
private double Determinant(double[,] a, double n)
{
int i, j, j1, j2;
double d = 0;
double[,] m =
{
{ ZERO, ZERO, ZERO, ZERO },
{ ZERO, ZERO, ZERO, ZERO },
{ ZERO, ZERO, ZERO, ZERO },
{ ZERO, ZERO, ZERO, ZERO }
};
if (n == 2)
{
// Terminate recursion.
d = a[0, 0] * a[1, 1] - a[1, 0] * a[0, 1];
}
else
{
d = 0;
for (j1 = 0; j1 < n; j1++) // Do each column.
{
for (i = 1; i < n; i++) // Create minor.
{
j2 = 0;
for (j = 0; j < n; j++)
{
if (j == j1) continue;
m[i - 1, j2] = a[i, j];
j2++;
}
}
// Sum (+/-)cofactor * minor.
d = d + System.Math.Pow(-1.0, j1) * a[0, j1] * this.Determinant(m, n - 1);
}
}
return d;
}
}
我是怎么说的，我的数据球体数量可能会有所不同，但我最多有 36 个球体，每个球体有 4 个点 x、y、z。 如果我可以将结果中心存储在另一个列表中，这将非常有用，也许类似于：
ballCent[0]=
center-> [0][1][2] //center of the sphere[x][y][z].
radius-> [0] //radius of the sphere.
我希望我解释得足够清楚，我不是英语母语人士，我非常感谢社区的帮助。 附言。我个人用我的数据尝试了该程序的java版本，它非常适合我。 链接在这里： http://www.convertalot.com/sphere_solver.html
@Robert Bruce 出色的 C++ 答案的 C# 端口。我并不声称自己理解数学。但我确实添加了一个测试用例并且它有效！
using System;
/*
https://stackoverflow.com/a/70846528/230851
class Point {
public:
double x;
double y;
double z;
Point() { x = 0; y = 0; z = 0; }
Point(double x_, double y_, double z_) { x = x_; y = y_; z = z_; }
};
class Sphere {
public:
Point center;
double radius;
Sphere(Point center_, double radius_) {
center = Point(center_.x, center_.y, center_.z);
radius = radius_;
}
};
Sphere
sphereFromFourPoints(Point a, Point b, Point c, Point d)
{
#define U(a,b,c,d,e,f,g,h) (a.z - b.z)*(c.x*d.y - d.x*c.y) - (e.z - f.z)*(g.x*h.y - h.x*g.y)
#define D(x,y,a,b,c) (a.x*(b.y-c.y) + b.x*(c.y-a.y) + c.x*(a.y-b.y))
#define E(x,y) ((ra*D(x,y,b,c,d) - rb*D(x,y,c,d,a) + rc*D(x,y,d,a,b) - rd*D(x,y,a,b,c)) / uvw)
double u = U(a,b,c,d,b,c,d,a);
double v = U(c,d,a,b,d,a,b,c);
double w = U(a,c,d,b,b,d,a,c);
double uvw = 2 * (u + v + w);
if (uvw == 0.0) {
// Oops. The points are coplanar.
}
auto sq = [] (Point p) { return p.x*p.x + p.y*p.y + p.z*p.z; };
double ra = sq(a);
double rb = sq(b);
double rc = sq(c);
double rd = sq(d);
double x0 = E(y,z);
double y0 = E(z,x);
double z0 = E(x,y);
double radius = sqrt(sq(Point(a.x - x0, a.y - y0, a.z - z0)));
return Sphere(Point(x0, y0, z0), radius);
}
*/
static class SphereSolver
{
public struct Point
{
public double x;
public double y;
public double z;
public Point(double x, double y, double z)
{
this.x = x;
this.y = y;
this.z = z;
}
public double component(int n)
{
switch (n)
{
case 0: return x;
case 1: return y;
case 2: return z;
default: throw new Exception();
}
}
};
public struct Sphere
{
public Point center;
public double radius;
public Sphere(Point center, double radius)
{
this.center = center;
this.radius = radius;
}
};
public static Sphere SphereFromFourPoints(Point a, Point b, Point c, Point d)
{
static double U(Point a, Point b, Point c, Point d, Point e, Point f, Point g, Point h)
{
return (a.z - b.z) * (c.x * d.y - d.x * c.y) - (e.z - f.z) * (g.x * h.y - h.x * g.y);
}
static double D(int x, int y, Point a, Point b, Point c)
{
return a.component(x) * (b.component(y) - c.component(y)) +
b.component(x) * (c.component(y) - a.component(y)) +
c.component(x) * (a.component(y) - b.component(y));
}
static double E(int x, int y, Point a, Point b, Point c, Point d, double ra, double rb, double rc, double rd, double uvw)
{
return ( ra * D(x, y, b, c, d) - rb * D(x, y, c, d, a) +
rc * D(x, y, d, a, b) - rd * D(x, y, a, b, c) ) / uvw;
}
double u = U(a, b, c, d, b, c, d, a);
double v = U(c, d, a, b, d, a, b, c);
double w = U(a, c, d, b, b, d, a, c);
double uvw = 2 * (u + v + w);
if (uvw == 0.0)
{
// Oops. The points are coplanar.
// You probably want to replace this with abs(uvw) < epsilon, with some epsilon appropriate for your project.
}
static double sq(Point p)
{
return p.x * p.x + p.y * p.y + p.z * p.z;
}
int x = 0;
int y = 1;
int z = 2;
double ra = sq(a);
double rb = sq(b);
double rc = sq(c);
double rd = sq(d);
double x0 = E(y, z, a, b, c, d, ra, rb, rc, rd, uvw);
double y0 = E(z, x, a, b, c, d, ra, rb, rc, rd, uvw);
double z0 = E(x, y, a, b, c, d, ra, rb, rc, rd, uvw);
double radius = System.Math.Sqrt(sq(new Point(a.x - x0, a.y - y0, a.z - z0)));
return new Sphere(new Point(x0, y0, z0), radius);
}
}
class Program
{
static void Main(string[] args)
{
SphereSolver.Sphere targetSphere = new SphereSolver.Sphere(new SphereSolver.Point(1, 2, 3), 5);
SphereSolver.Point[] randomDirections =
{
new SphereSolver.Point(1, 1, 1),
new SphereSolver.Point(7, 1, 5),
new SphereSolver.Point(0, -4, -3),
new SphereSolver.Point(-3, -5, -2),
};
SphereSolver.Point[] pointsOnTestSphere = new SphereSolver.Point[4];
for (int n = 0; n < pointsOnTestSphere.Length; ++n)
{
SphereSolver.Point d = randomDirections[n];
SphereSolver.Point p = new SphereSolver.Point();
double length = Math.Sqrt(d.x * d.x + d.y * d.y + d.z * d.z);
p.x = d.x * targetSphere.radius / length;
p.y = d.y * targetSphere.radius / length;
p.z = d.z * targetSphere.radius / length;
p.x += targetSphere.center.x;
p.y += targetSphere.center.y;
p.z += targetSphere.center.z;
pointsOnTestSphere[n] = p;
}
SphereSolver.Sphere answerSphere = SphereSolver.SphereFromFourPoints(
pointsOnTestSphere[0],
pointsOnTestSphere[1],
pointsOnTestSphere[2],
pointsOnTestSphere[3]
);
Console.WriteLine("Target Sphere: {0}, {1}, {2}, radius {3}", targetSphere.center.x, targetSphere.center.y, targetSphere.center.z, targetSphere.radius);
Console.WriteLine("Answer Sphere: {0}, {1}, {2}, radius {3}", answerSphere.center.x, answerSphere.center.y, answerSphere.center.z, answerSphere.radius);
}
}
给定四个点，a、b、c和d，您可以通过将以下行列式设置为零并求解来找到中心：
| (x^2 + y^2 + z^2) x y z 1 |
| (ax^2 + ay^2 + az^2) ax ay az 1 |
| (bx^2 + by^2 + bz^2) bx by bz 1 | = 0.
| (cx^2 + cy^2 + cz^2) cx cy cz 1 |
| (dx^2 + dy^2 + dz^2) dx dy dz 1 |
数学很粗糙，但以下 C++ 代码实现了解决方案：
class Point {
public:
double x;
double y;
double z;
Point() { x = 0; y = 0; z = 0; }
Point(double x_, double y_, double z_) { x = x_; y = y_; z = z_; }
};
class Sphere {
public:
Point center;
double radius;
Sphere(Point center_, double radius_) {
center = Point(center_.x, center_.y, center_.z);
radius = radius_;
}
};
Sphere
sphereFromFourPoints(Point a, Point b, Point c, Point d)
{
#define U(a,b,c,d,e,f,g,h) (a.z - b.z)*(c.x*d.y - d.x*c.y) - (e.z - f.z)*(g.x*h.y - h.x*g.y)
#define D(x,y,a,b,c) (a.x*(b.y-c.y) + b.x*(c.y-a.y) + c.x*(a.y-b.y))
#define E(x,y) ((ra*D(x,y,b,c,d) - rb*D(x,y,c,d,a) + rc*D(x,y,d,a,b) - rd*D(x,y,a,b,c)) / uvw)
double u = U(a,b,c,d,b,c,d,a);
double v = U(c,d,a,b,d,a,b,c);
double w = U(a,c,d,b,b,d,a,c);
double uvw = 2 * (u + v + w);
if (uvw == 0.0) {
// Oops. The points are coplanar.
}
auto sq = [] (Point p) { return p.x*p.x + p.y*p.y + p.z*p.z; };
double ra = sq(a);
double rb = sq(b);
double rc = sq(c);
double rd = sq(d);
double x0 = E(y,z);
double y0 = E(z,x);
double z0 = E(x,y);
double radius = sqrt(sq(Point(a.x - x0, a.y - y0, a.z - z0)));
return Sphere(Point(x0, y0, z0), radius);
#undef U
#undef D
#undef E
}
我编写了这段代码，并特此将其置于公共领域。你可以用它做任何你想做的事。
这是一个很容易理解的几何解：
注意：有一些众所周知的方法可以从 3 个点得到圆。 圆的“轴”是穿过圆心且垂直于圆平面的线。
使用前两个点定义一条线。对后两个做同样的事情。找到两条线之间的 POI。