From 39744e9ea26132fa78a98017aa8a8add22aededb Mon Sep 17 00:00:00 2001 From: sangminpakr <82554467+sangminpakr@users.noreply.github.com> Date: Tue, 27 Apr 2021 19:09:08 +0900 Subject: [PATCH] Add sphere3D --- src/Spatial.Tests/Euclidean/Sphere3DTests.cs | 84 +++++++ src/Spatial/Euclidean/Sphere3D.cs | 245 +++++++++++++++++++ 2 files changed, 329 insertions(+) create mode 100644 src/Spatial.Tests/Euclidean/Sphere3DTests.cs create mode 100644 src/Spatial/Euclidean/Sphere3D.cs diff --git a/src/Spatial.Tests/Euclidean/Sphere3DTests.cs b/src/Spatial.Tests/Euclidean/Sphere3DTests.cs new file mode 100644 index 00000000..f1cdd0b9 --- /dev/null +++ b/src/Spatial.Tests/Euclidean/Sphere3DTests.cs @@ -0,0 +1,84 @@ + +using System; +using MathNet.Spatial.Euclidean; +using NUnit.Framework; + +namespace MathNet.Spatial.UnitTests.Euclidean +{ + [TestFixture] + public class Sphere3DTests + { + [TestCase("2.3,4,5", 1)] + [TestCase("-1,4,2", 1)] + [TestCase("3,-2,1", 1)] + [TestCase("2,4,-4", 1)] + [TestCase("0,0,0", 1)] + [TestCase("0,0,0", 2)] + public void SphereCenterRadius(string p1s, double radius) + { + var center = Point3D.Parse(p1s); + var sphere = new Sphere3D(center, radius); + Assert.AreEqual(2 * radius, sphere.Diameter, double.Epsilon); + Assert.AreEqual(2 * Math.PI * radius, sphere.Circumference, double.Epsilon); + Assert.AreEqual(Math.PI * radius * radius * 4, sphere.Area, double.Epsilon); + } + + [TestCase("-1,0,0", "1,0,0","0,0,0", 1)] + [TestCase("-1,0,0", "5,0,0","2,0,0", 3)] + [TestCase("0,0,0", "0,2,0", "0,1,0", 1)] + [TestCase("2,4,0", "2,2,2", "2,3,1", 1.4142135623)] + [TestCase("1,2,0", "0,0,0", "0.5,1,0", 1.11803398874)] + [TestCase("2,3,4", "4,5,2", "3,4,3", 1.73205080)] + public void SphereFromTwoPoints(string p1s, string p2s, string p2centers, double p2radius) + { + var p2p1 = Point3D.Parse(p1s); + var p2p2 = Point3D.Parse(p2s); + var p2center = Point3D.Parse(p2centers); + var p2sphere = Sphere3D.FromTwoPoints(p2p1, p2p2); + + AssertGeometry.AreEqual(p2center, p2sphere.CenterPoint); + Assert.AreEqual(p2radius, p2sphere.Radius, 1e-6); + } + + + [TestCase("3,2,1", "1,-2,-3", "2,1,3", "-1,1,2", "1.2631578, -0.8421052, 0.2105263", 3.4230763)] + [TestCase("1,2,1", "0,3,2", "1,1,0", "-1,0,-2", "-1.5,2.5,-0.5", 2.9580398)] + [TestCase("2,1,1", "5,2,2", "2,3,3", "3,6,2", "3.1666666,3.5833333,0.4166666", 2.8939592)] + [TestCase("2,1,1", "4,2,2", "2,1,3", "1,3,2", "2.4,2.2,2", 1.612451549)] + [TestCase("4,2,1", "2,-2,-3", "6,1,3", "-1,1,2", "2.5405405,-2.986486486,2.21621621621", 5.336126992333)] + public void SphereFromFourPoints(string p1s, string p2s, string p3s, string p4s, string centers, double radius) + { + var p1 = Point3D.Parse(p1s); + var p2 = Point3D.Parse(p2s); + var p3 = Point3D.Parse(p3s); + var p4 = Point3D.Parse(p4s); + var center = Point3D.Parse(centers); + var sphere = Sphere3D.FromFourPoints(p1, p2, p3, p4); + + AssertGeometry.AreEqual(center, sphere.CenterPoint); + Assert.AreEqual(radius, sphere.Radius, 1e-6); + } + + [Test] + public void SphereFromFourPointsArgumentException() + { + var p1 = new Point3D(0, 0, 0); + var p2 = new Point3D(1, 2, 4); + var p3 = new Point3D(2, 4, 8); + var p4 = new Point3D(-1, 1, 2); + + Assert.Throws(() => Sphere3D.FromFourPoints(p1, p2, p3, p4)); + } + + [Test] + public void SphereFromTwoPointsArgumentException() + { + var p1 = new Point3D(1, 4, 1); + var p2 = new Point3D(1, 4, 1); + + Assert.Throws(() => Sphere3D.FromTwoPoints(p1, p2)); + } + } +} + + diff --git a/src/Spatial/Euclidean/Sphere3D.cs b/src/Spatial/Euclidean/Sphere3D.cs new file mode 100644 index 00000000..f94727d2 --- /dev/null +++ b/src/Spatial/Euclidean/Sphere3D.cs @@ -0,0 +1,245 @@ +using System; +using System.Diagnostics.Contracts; +using MathNet.Numerics.LinearAlgebra; +using HashCode = MathNet.Spatial.Internals.HashCode; + +namespace MathNet.Spatial.Euclidean +{ + /// + /// Describes a 3 dimensional sphere. + /// + [Serializable] + public struct Sphere3D : IEquatable + { + /// + /// The center of the sphere. + /// + public readonly Point3D CenterPoint; + + /// + /// The radius of the sphere. + /// + public readonly double Radius; + + /// + /// Initializes a new instance of the . + /// + /// The center of the Sphere + /// The radius of the Sphere + public Sphere3D(Point3D centerPoint, double radius) + { + // Anypoint is same to other point. + if (radius <= 0) + { + throw new ArgumentException("The radius is negative."); + } + + this.CenterPoint = centerPoint; + this.Radius = radius; + } + + /// + /// Gets the diameter of the sphere. + /// + [Pure] + public double Diameter => 2 * Radius; + + /// + /// Gets the circumference of the sphere. + /// + [Pure] + public double Circumference => 2 * Math.PI * Radius; + + /// + /// Gets the volume of the sphere. + /// + [Pure] + public double Volume => 4 / 3 * Math.PI * this.Radius * this.Radius * this.Radius; + + /// + /// Gets the area of the sphere. + /// + [Pure] + public double Area => 4 * Math.PI * this.Radius * this.Radius; + + /// + /// Returns a value that indicates whether each pair of elements in two specified sphere is equal. + /// + /// The first sphere to compare + /// The second sphere to compare + /// True if the sphere are the same; otherwise false. + public static bool operator ==(Sphere3D left, Sphere3D right) + { + return left.Equals(right); + } + + /// + /// Returns a value that indicates whether any pair of elements in two specified sphere is not equal. + /// + /// The first sphere to compare + /// The second sphere to compare + /// True if the sphere are different; otherwise false. + public static bool operator !=(Sphere3D left, Sphere3D right) + { + return !left.Equals(right); + } + + /// + /// Creates an instance of sphere which lie along its circumference of tetrahydron made of four points. + /// + /// The first point on the sphere + /// The second point on the sphere + /// The third point on the sphere + /// The last point on the sphere + public static Sphere3D FromFourPoints(Point3D p1, Point3D p2, Point3D p3, Point3D p4) + { + // https://mathworld.wolfram.com/Circumsphere.html + // + // The equation for the circumsphere of the tetrahedron with polygon vertices p1 = { x1, y1, z1 }, p2 = { x2, y2, z2 }, p3 = { x3, y3, z3 }, p4 = { x4, y4, z4 }, + // | x1^2 + y1^2 + z1^2 x1 y1 z1 1 | = 0 + // | x2^2 + y2^2 + z2^2 x2 y2 z2 1 | + // | x3^2 + y3^2 + z3^2 x3 y3 z3 1 | + // | x4^2 + y4^2 + z4^2 x4 y4 z4 1 | + // + // gives + // a(x^2 + y^2 + z^2) - (Dx*x + Dy*y + Dz*z) + c = 0, + // + // where + // A = | x1 y2 z3 1 | = 0 + // | x2 y2 z2 1 | + // | x3 y3 z3 1 | + // | x4 y4 z4 1 |, + // and + // Dx = | x1^2+y1^2+z1^2 y1 z1 1 | Dy= - | x1^2+y1^2+z1^2x1z1 1 | Dz = | x1^2+y1^2+z1^2 x1 y1 1 | c= | x1^2+y1^2+z1^2 x1 y1 z1 | + // | x2^2+y2^2+z2^2 y2 z2 1 | | x2^2+y2^2+z2^2x2z2 1 | | x2^2+y2^2+z2^2 x2 y2 1 | | x2^2+y2^2+z2^2 x2 y2 z2 | + // | x3^2+y3^2+z3^2 y3 z3 1 | | x3^2+y3^2+z3^2x3z3 1 | | x3^2+y3^2+z3^2 x3 y3 1 | | x3^2+y3^2+z3^2 x3 y3 z3 | + // | x4^2+y4^2+z4^2 y4 z4 1 | | x4^2+y4^2+z4^2x4z4 1 | | x4^2+y4^2+z4^2 x4 y4 1 | | x4^2+y4^2+z4^2 x4 y4 z4 | . + // + // Completing the square gives, + // a(x-(Dx)/(2a))^2 + a(y-(Dy)/(2a))^2 + a(z-(Dz)/(2a))^2 - (Dx^2+Dy^2+Dz^2)/(4a)+c = 0, + // which is a sphere of the form, + // (x-x0)^2 + (y-y0)^2 + (z-z0)^2 = r^2, + // + // With circumcenter and circumradius. + // x0 = (Dx)/(2a), y0 = (Dy)/(2a), z0 =(Dz)/(2a), and r = (sqrt(Dx^2+Dy^2+Dz^2-4ac))/(2|a|). + + if (p1 == p2 || p1 == p3 || p1 == p4 || p2 == p3 || p2 == p4 || p3 ==p4) + { + throw new ArgumentException("All points must be differrent."); + } + + var Dx = Matrix.Build.DenseOfArray(new[,] + { + { p1.X*p1.X + p1.Y*p1.Y + p1.Z*p1.Z, p1.Y,p1.Z, 1 }, + { p2.X*p2.X + p2.Y*p2.Y + p2.Z*p2.Z, p2.Y,p2.Z, 1 }, + { p3.X*p3.X + p3.Y*p3.Y + p3.Z*p3.Z, p3.Y,p3.Z, 1 }, + { p4.X*p4.X + p4.Y*p4.Y + p4.Z*p4.Z, p4.Y,p4.Z, 1 } + }); + var Dy = Matrix.Build.DenseOfArray(new[,] + { + { p1.X*p1.X + p1.Y*p1.Y + p1.Z*p1.Z, p1.X,p1.Z, 1 }, + { p2.X*p2.X + p2.Y*p2.Y + p2.Z*p2.Z, p2.X,p2.Z, 1 }, + { p3.X*p3.X + p3.Y*p3.Y + p3.Z*p3.Z, p3.X,p3.Z, 1 }, + { p4.X*p4.X + p4.Y*p4.Y + p4.Z*p4.Z, p4.X,p4.Z, 1 } + }); + var Dz = Matrix.Build.DenseOfArray(new[,] + { + { p1.X*p1.X + p1.Y*p1.Y + p1.Z*p1.Z, p1.X,p1.Y, 1 }, + { p2.X*p2.X + p2.Y*p2.Y + p2.Z*p2.Z, p2.X,p2.Y ,1 }, + { p3.X*p3.X + p3.Y*p3.Y + p3.Z*p3.Z, p3.X,p3.Y, 1 }, + { p4.X*p4.X + p4.Y*p4.Y + p4.Z*p4.Z, p4.X,p4.Y, 1 } + }); + var a = Matrix.Build.DenseOfArray(new[,] + { + { p1.X, p1.Y, p1.Z, 1 }, + { p2.X, p2.Y, p2.Z, 1 }, + { p3.X, p3.Y, p3.Z, 1 }, + { p4.X, p4.Y, p4.Z, 1 } + }); + var c = Matrix.Build.DenseOfArray(new[,] + { + { p1.X*p1.X + p1.Y*p1.Y + p1.Z*p1.Z, p1.X, p1.Y, p1.Z }, + { p2.X*p2.X + p2.Y*p2.Y + p2.Z*p2.Z, p2.X, p2.Y, p2.Z }, + { p3.X*p3.X + p3.Y*p3.Y + p3.Z*p3.Z, p3.X, p3.Y, p3.Z }, + { p4.X*p4.X + p4.Y*p4.Y + p4.Z*p4.Z, p4.X, p4.Y, p4.Z } + }); + + var detDx = Dx.Determinant(); + var detDy = -Dy.Determinant(); + var detDz = Dz.Determinant(); + var deta = a.Determinant(); + var detc = c.Determinant(); + + if (deta == 0) + { + throw new ArgumentException("A circumcenter cannot be created from these points, are they collinear?"); + } + + // Finally, we can get the circumcenter and radius. + var x0 = detDx / (2 * deta); + var y0 = detDy / (2 * deta); + var z0 = detDz / (2 * deta); + var center = new Point3D(x0, y0, z0); + var radius = Math.Sqrt(detDx * detDx + detDy * detDy + detDz * detDz - 4 * deta * detc) / (2 * Math.Abs(deta)); + + return new Sphere3D(center,radius); + } + + /// + /// Creates an instance of from 2 points of diameter. + /// + /// The start point of diameter. + /// The end point of diameter. + /// + public static Sphere3D FromTwoPoints(Point3D p1, Point3D p2) + { + if (p1 == p2) + { + throw new ArgumentException("A sphere cannot be created from two identical points."); + } + + var center = Point3D.MidPoint(p1, p2); + var radius = center.DistanceTo(p1); + + return new Sphere3D(center, radius); + } + + /// + /// Returns a value to indicate if a pair of sphere are equal. + /// + /// The sphere to compare against. + /// A tolerance (epsilon) to adjust for floating point error. + /// True if the points are equal; otherwise false. + public bool Equals(Sphere3D c, double tolerance) + { + if (tolerance < 0) + { + throw new ArgumentException("epsilon < 0"); + } + + return Math.Abs(c.Radius - this.Radius) < tolerance + && this.CenterPoint.Equals(c.CenterPoint, tolerance); + } + + /// + [Pure] + public bool Equals(Sphere3D c) + { + return this.CenterPoint.Equals(c.CenterPoint) + && this.Radius.Equals(c.Radius); + } + + /// + [Pure] + public override bool Equals(object obj) => obj is Sphere3D c && this.Equals(c); + + /// + [Pure] + public override int GetHashCode() => HashCode.Combine(this.CenterPoint, this.Radius); + } + + +} + + +