如何用4个点计算球心?

问题描述 投票:0回答:4

问题是我对编程非常陌生,现在我需要制作一个程序来计算许多球体的中心(最多 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

c# class math determinants
4个回答
1
投票

@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);
    }
}

1
投票

给定四个点,abcd,您可以通过将以下行列式设置为零并求解来找到中心:

| (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
}

我编写了这段代码,并特此将其置于公共领域。你可以用它做任何你想做的事。


0
投票

这是一个很容易理解的几何解:

  1. 选择 3 个点并获得穿过它们的唯一 3D 圆 C1。
  2. 选择另外 3 个点并获得另一个穿过它们的圆 C2。
  3. 获取C1的轴线L1,C2的轴线L2
  4. 找到 L1 和 L2 在 3D 中交叉的位置

注意:有一些众所周知的方法可以从 3 个点得到圆。 圆的“轴”是穿过圆心且垂直于圆平面的线。


-1
投票

使用前两个点定义一条线。对后两个做同样的事情。找到两条线之间的 POI。

© www.soinside.com 2019 - 2024. All rights reserved.