From 3611b34a4ce5a4f9d1072711b1e1a7e611e91483 Mon Sep 17 00:00:00 2001 From: allem Date: Tue, 3 Apr 2018 23:12:35 -0400 Subject: [PATCH 1/7] mvec4d with inverse --- Tetrahedrons/MVec3d.cs | 340 ++++++----- Tetrahedrons/MVec4D.cs | 970 ++++++++++++++++++++++++++++++ Tetrahedrons/Program.cs | 33 +- Tetrahedrons/TetrahedronWindow.cs | 38 +- Tetrahedrons/Tetrahedrons.csproj | 1 + 5 files changed, 1183 insertions(+), 199 deletions(-) create mode 100644 Tetrahedrons/MVec4D.cs diff --git a/Tetrahedrons/MVec3d.cs b/Tetrahedrons/MVec3d.cs index 97bf346..4ce9eb8 100644 --- a/Tetrahedrons/MVec3d.cs +++ b/Tetrahedrons/MVec3d.cs @@ -3,205 +3,203 @@ using OpenTK; namespace Tetrahedrons { - // ReSharper disable once InconsistentNaming - public struct MVec3d - { - public static readonly MVec3d Zero = new MVec3d(0, 0, 0, 0, 0, 0, 0, 0); - public static readonly MVec3d One = new MVec3d(1, 1, 1, 1, 1, 1, 1, 1); + public struct MVec3D + { + public static readonly MVec3D Zero = new MVec3D(0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec3D One = new MVec3D(1, 1, 1, 1, 1, 1, 1, 1); - public static readonly MVec3d Unit = new MVec3d(1, 0, 0, 0, 0, 0, 0, 0); - public static readonly MVec3d Unit1 = new MVec3d(0, 1, 0, 0, 0, 0, 0, 0); - public static readonly MVec3d Unit2 = new MVec3d(0, 0, 1, 0, 0, 0, 0, 0); - public static readonly MVec3d Unit3 = new MVec3d(0, 0, 0, 1, 0, 0, 0, 0); - public static readonly MVec3d Unit23 = new MVec3d(0, 0, 0, 0, 1, 0, 0, 0); - public static readonly MVec3d Unit31 = new MVec3d(0, 0, 0, 0, 0, 1, 0, 0); - public static readonly MVec3d Unit12 = new MVec3d(0, 0, 0, 0, 0, 0, 1, 0); - public static readonly MVec3d Unit123 = new MVec3d(0, 0, 0, 0, 0, 0, 0, 1); + public static readonly MVec3D Unit = new MVec3D(1, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec3D Unit1 = new MVec3D(0, 1, 0, 0, 0, 0, 0, 0); + public static readonly MVec3D Unit2 = new MVec3D(0, 0, 1, 0, 0, 0, 0, 0); + public static readonly MVec3D Unit3 = new MVec3D(0, 0, 0, 1, 0, 0, 0, 0); + public static readonly MVec3D Unit23 = new MVec3D(0, 0, 0, 0, 1, 0, 0, 0); + public static readonly MVec3D Unit31 = new MVec3D(0, 0, 0, 0, 0, 1, 0, 0); + public static readonly MVec3D Unit12 = new MVec3D(0, 0, 0, 0, 0, 0, 1, 0); + public static readonly MVec3D Unit123 = new MVec3D(0, 0, 0, 0, 0, 0, 0, 1); - public double E; - public double E1; - public double E2; - public double E3; - public double E23; - public double E31; - public double E12; - public double E123; + public double E; + public double E1; + public double E2; + public double E3; + public double E23; + public double E31; + public double E12; + public double E123; - public MVec3d(double e, double e1, double e2, double e3, double e23, double e31, double e12, double e123) - { - E = e; - E1 = e1; - E2 = e2; - E3 = e3; - E23 = e23; - E31 = e31; - E12 = e12; - E123 = e123; - } + public MVec3D(double e, double e1, double e2, double e3, double e23, double e31, double e12, double e123) + { + E = e; + E1 = e1; + E2 = e2; + E3 = e3; + E23 = e23; + E31 = e31; + E12 = e12; + E123 = e123; + } - public MVec3d(MVec3d s, MVec3d v, MVec3d b, MVec3d t) - : this(s.E, v.E1, v.E2, v.E3, b.E23, b.E31, b.E12, t.E123) - { - } + public MVec3D(MVec3D s, MVec3D v, MVec3D b, MVec3D t) + : this(s.E, v.E1, v.E2, v.E3, b.E23, b.E31, b.E12, t.E123) + { + } - public static MVec3d Scalar(double e) - { - return new MVec3d(e, 0, 0, 0, 0, 0, 0, 0); - } + public static MVec3D Scalar(double e) + { + return new MVec3D(e, 0, 0, 0, 0, 0, 0, 0); + } - public static MVec3d Vector(double e1, double e2, double e3) - { - return new MVec3d(0, e1, e2, e3, 0, 0, 0, 0); - } + public static MVec3D Vector(double e1, double e2, double e3) + { + return new MVec3D(0, e1, e2, e3, 0, 0, 0, 0); + } - public static MVec3d Vector(Vector3d v) - { - return new MVec3d(0, v.X, v.Y, v.Z, 0, 0, 0, 0); - } + public static MVec3D Vector(Vector3d v) + { + return new MVec3D(0, v.X, v.Y, v.Z, 0, 0, 0, 0); + } - public static MVec3d Bivector(double e23, double e31, double e12) - { - return new MVec3d(0, 0, 0, 0, e23, e31, e12, 0); - } + public static MVec3D Bivector(double e23, double e31, double e12) + { + return new MVec3D(0, 0, 0, 0, e23, e31, e12, 0); + } - public static MVec3d Trivector(double e123) - { - return new MVec3d(0, 0, 0, 0, 0, 0, 0, e123); - } + public static MVec3D Trivector(double e123) + { + return new MVec3D(0, 0, 0, 0, 0, 0, 0, e123); + } - public static MVec3d Complex(double real, double imag) - { - return new MVec3d(real, 0, 0, 0, 0, 0, 0, imag); - } + public static MVec3D Complex(double real, double imag) + { + return new MVec3D(real, 0, 0, 0, 0, 0, 0, imag); + } - public static MVec3d Complex(double w, double i, double j, double k) - { - return new MVec3d(w, 0, 0, 0, i, j, k, 0); - } + public static MVec3D Complex(double w, double i, double j, double k) + { + return new MVec3D(w, 0, 0, 0, i, j, k, 0); + } - public static MVec3d Complex(double angle, MVec3d plane) - { - return Math.Cos(angle) + Math.Sin(angle) * plane; - } + public static MVec3D Complex(double angle, MVec3D plane) + { + return Math.Cos(angle) + Math.Sin(angle) * plane; + } - public static MVec3d Add(MVec3d m, MVec3d n) - { - return new MVec3d(m.E + n.E, - m.E1 + n.E1, m.E2 + n.E2, m.E3 + n.E3, - m.E23 + n.E23, m.E31 + n.E31, m.E12 + n.E12, - m.E123 + n.E123); - } + public static MVec3D Add(MVec3D m, MVec3D n) + { + return new MVec3D(m.E + n.E, + m.E1 + n.E1, m.E2 + n.E2, m.E3 + n.E3, + m.E23 + n.E23, m.E31 + n.E31, m.E12 + n.E12, + m.E123 + n.E123); + } - public static MVec3d Sub(MVec3d m, MVec3d n) - { - return new MVec3d(m.E - n.E, - m.E1 - n.E1, m.E2 - n.E2, m.E3 - n.E3, - m.E23 - n.E23, m.E31 - n.E31, m.E12 - n.E12, - m.E123 - n.E123); - } + public static MVec3D Sub(MVec3D m, MVec3D n) + { + return new MVec3D(m.E - n.E, + m.E1 - n.E1, m.E2 - n.E2, m.E3 - n.E3, + m.E23 - n.E23, m.E31 - n.E31, m.E12 - n.E12, + m.E123 - n.E123); + } - public static MVec3d Mult(MVec3d m, MVec3d n) - { - return new MVec3d( - m.E * n.E + m.E1 * n.E1 + m.E2 * n.E2 + m.E3 * n.E3 - m.E23 * n.E23 - m.E31 * n.E31 - m.E12 * n.E12 - m.E123 * n.E123, - m.E * n.E1 + m.E1 * n.E - m.E2 * n.E12 + m.E3 * n.E31 - m.E23 * n.E123 - m.E31 * n.E3 + m.E12 * n.E2 - m.E123 * n.E23, - m.E * n.E2 + m.E1 * n.E12 + m.E2 * n.E - m.E3 * n.E23 + m.E23 * n.E3 - m.E31 * n.E123 - m.E12 * n.E1 - m.E123 * n.E31, - m.E * n.E3 - m.E1 * n.E31 + m.E2 * n.E23 + m.E3 * n.E - m.E23 * n.E2 + m.E31 * n.E1 - m.E12 * n.E123 - m.E123 * n.E12, - m.E * n.E23 + m.E1 * n.E123 + m.E2 * n.E3 - m.E3 * n.E2 + m.E23 * n.E - m.E31 * n.E12 + m.E12 * n.E31 + m.E123 * n.E1, - m.E * n.E31 - m.E1 * n.E3 + m.E2 * n.E123 + m.E3 * n.E1 + m.E23 * n.E12 + m.E31 * n.E - m.E12 * n.E23 + m.E123 * n.E2, - m.E * n.E12 + m.E1 * n.E2 - m.E2 * n.E1 + m.E3 * n.E123 - m.E23 * n.E31 + m.E31 * n.E23 + m.E12 * n.E + m.E123 * n.E3, - m.E * n.E123 + m.E1 * n.E23 + m.E2 * n.E31 + m.E3 * n.E12 + m.E23 * n.E1 + m.E31 * n.E2 + m.E12 * n.E3 + m.E123 * n.E - ); - } + public static MVec3D Mult(MVec3D m, MVec3D n) + { + return new MVec3D( + m.E * n.E + m.E1 * n.E1 + m.E2 * n.E2 + m.E3 * n.E3 - m.E23 * n.E23 - m.E31 * n.E31 - m.E12 * n.E12 - m.E123 * n.E123, + m.E * n.E1 + m.E1 * n.E - m.E2 * n.E12 + m.E3 * n.E31 - m.E23 * n.E123 - m.E31 * n.E3 + m.E12 * n.E2 - m.E123 * n.E23, + m.E * n.E2 + m.E1 * n.E12 + m.E2 * n.E - m.E3 * n.E23 + m.E23 * n.E3 - m.E31 * n.E123 - m.E12 * n.E1 - m.E123 * n.E31, + m.E * n.E3 - m.E1 * n.E31 + m.E2 * n.E23 + m.E3 * n.E - m.E23 * n.E2 + m.E31 * n.E1 - m.E12 * n.E123 - m.E123 * n.E12, + m.E * n.E23 + m.E1 * n.E123 + m.E2 * n.E3 - m.E3 * n.E2 + m.E23 * n.E - m.E31 * n.E12 + m.E12 * n.E31 + m.E123 * n.E1, + m.E * n.E31 - m.E1 * n.E3 + m.E2 * n.E123 + m.E3 * n.E1 + m.E23 * n.E12 + m.E31 * n.E - m.E12 * n.E23 + m.E123 * n.E2, + m.E * n.E12 + m.E1 * n.E2 - m.E2 * n.E1 + m.E3 * n.E123 - m.E23 * n.E31 + m.E31 * n.E23 + m.E12 * n.E + m.E123 * n.E3, + m.E * n.E123 + m.E1 * n.E23 + m.E2 * n.E31 + m.E3 * n.E12 + m.E23 * n.E1 + m.E31 * n.E2 + m.E12 * n.E3 + m.E123 * n.E + ); + } - public double Norm2 => E * E + E1 * E1 + E2 * E2 + E3 * E3 + E23 * E23 + E31 * E31 + E12 * E12 + E123 * E123; - public double Norm => Math.Sqrt(Norm2); + public double Norm2 => E * E + E1 * E1 + E2 * E2 + E3 * E3 + E23 * E23 + E31 * E31 + E12 * E12 + E123 * E123; + public double Norm => Math.Sqrt(Norm2); - public MVec3d Inv - { - get - { - var n2 = Norm2; - return new MVec3d(E / n2, -E1 / n2, -E2 / n2, -E3 / n2, -E23 / n2, -E31 / n2, -E12 / n2, E123 / n2); -// return new MVec3d(E / n2, E1 / n2, E2 / n2, E3 / n2, E23 / n2, E31 / n2, E12 / n2, E123 / n2); - } - } + public MVec3D Inv + { + get + { + var n2 = Norm2; + return new MVec3D(E / n2, -E1 / n2, -E2 / n2, -E3 / n2, -E23 / n2, -E31 / n2, -E12 / n2, E123 / n2); + } + } - public MVec3d Grade(int k) - { - switch (k) - { - case 0: - return Scalar(E); - case 1: - return Vector(E1, E2, E3); - case 2: - return Bivector(E23, E31, E12); - case 3: - return Trivector(E123); - default: - return new MVec3d(); - } - } + public MVec3D Grade(int k) + { + switch (k) + { + case 0: + return Scalar(E); + case 1: + return Vector(E1, E2, E3); + case 2: + return Bivector(E23, E31, E12); + case 3: + return Trivector(E123); + default: + return new MVec3D(); + } + } - public static MVec3d Inner(MVec3d m, MVec3d n) - { - var m0 = m.Grade(0); - var m1 = m.Grade(1); - var m2 = m.Grade(2); - var m3 = m.Grade(3); + public static MVec3D Inner(MVec3D m, MVec3D n) + { + var m0 = m.Grade(0); + var m1 = m.Grade(1); + var m2 = m.Grade(2); + var m3 = m.Grade(3); - var n0 = n.Grade(0); - var n1 = n.Grade(1); - var n2 = n.Grade(2); - var n3 = n.Grade(3); + var n0 = n.Grade(0); + var n1 = n.Grade(1); + var n2 = n.Grade(2); + var n3 = n.Grade(3); - return new MVec3d( - m0 * n0 + m1 * n1 + m2 * n2 + m3 * n3, - m0 * n1 + m1 * n2 + m2 * n3, - m0 * n2 + m1 * n3, - m0 * n3 - ); - } + return new MVec3D( + m0 * n0 + m1 * n1 + m2 * n2 + m3 * n3, + m0 * n1 + m1 * n2 + m2 * n3, + m0 * n2 + m1 * n3, + m0 * n3 + ); + } - public static MVec3d Outer(MVec3d m, MVec3d n) - { - var m0 = m.Grade(0); - var m1 = m.Grade(1); - var m2 = m.Grade(2); - var m3 = m.Grade(3); + public static MVec3D Outer(MVec3D m, MVec3D n) + { + var m0 = m.Grade(0); + var m1 = m.Grade(1); + var m2 = m.Grade(2); + var m3 = m.Grade(3); - var n0 = n.Grade(0); - var n1 = n.Grade(1); - var n2 = n.Grade(2); - var n3 = n.Grade(3); + var n0 = n.Grade(0); + var n1 = n.Grade(1); + var n2 = n.Grade(2); + var n3 = n.Grade(3); - return new MVec3d( - m0 * n0, - m1 * n0 + m0 * n1, - m2 * n0 + m1 * n1 + m0 * n2, - m3 * n0 + m2 * n1 + m1 * n2 + m0 * m3 - ); - } + return new MVec3D( + m0 * n0, + m1 * n0 + m0 * n1, + m2 * n0 + m1 * n1 + m0 * n2, + m3 * n0 + m2 * n1 + m1 * n2 + m0 * m3 + ); + } - public static MVec3d operator ~(MVec3d m) => m.Inv; - public static MVec3d operator *(MVec3d m, MVec3d n) => Mult(m, n); + public static MVec3D operator ~(MVec3D m) => m.Inv; + public static MVec3D operator *(MVec3D m, MVec3D n) => Mult(m, n); - public static MVec3d operator +(MVec3d m, MVec3d n) => Add(m, n); - public static MVec3d operator -(MVec3d m, MVec3d n) => Sub(m, n); + public static MVec3D operator +(MVec3D m, MVec3D n) => Add(m, n); + public static MVec3D operator -(MVec3D m, MVec3D n) => Sub(m, n); - public static MVec3d operator &(MVec3d m, MVec3d n) => Inner(m, n); - public static MVec3d operator ^(MVec3d m, MVec3d n) => Outer(m, n); + public static MVec3D operator &(MVec3D m, MVec3D n) => Inner(m, n); + public static MVec3D operator ^(MVec3D m, MVec3D n) => Outer(m, n); - public Vector3d VectorPart => new Vector3d(E1, E2, E3); + public Vector3d VectorPart => new Vector3d(E1, E2, E3); - public static implicit operator MVec3d(Vector3 v) => new MVec3d(0, v.X, v.Y, v.Z, 0, 0, 0, 0); - public static implicit operator MVec3d(Vector3d v) => new MVec3d(0, v.X, v.Y, v.Z, 0, 0, 0, 0); - public static implicit operator MVec3d(double s) => new MVec3d(s, 0, 0, 0, 0, 0, 0, 0); + public static implicit operator MVec3D(Vector3 v) => new MVec3D(0, v.X, v.Y, v.Z, 0, 0, 0, 0); + public static implicit operator MVec3D(Vector3d v) => new MVec3D(0, v.X, v.Y, v.Z, 0, 0, 0, 0); + public static implicit operator MVec3D(double s) => new MVec3D(s, 0, 0, 0, 0, 0, 0, 0); - public override string ToString() - { - return $"({E:0.00} {E1:0.00}e1 {E2:0.00}e2 {E3:0.00}e3 {E23:0.00}e23 {E31:0.00}e31 {E12:0.00}e12 {E123:0.00}e123)"; - } - } + public override string ToString() + { + return $"({E:0.00} {E1:0.00}e1 {E2:0.00}e2 {E3:0.00}e3 {E23:0.00}e23 {E31:0.00}e31 {E12:0.00}e12 {E123:0.00}e123)"; + } + } } \ No newline at end of file diff --git a/Tetrahedrons/MVec4D.cs b/Tetrahedrons/MVec4D.cs new file mode 100644 index 0000000..62216a0 --- /dev/null +++ b/Tetrahedrons/MVec4D.cs @@ -0,0 +1,970 @@ +using System; +using System.Collections.Generic; +using System.Runtime.Remoting.Messaging; +using System.Text; +using OpenTK; + +namespace Tetrahedrons +{ + public struct MVec4D + { + #region Units + + public static readonly MVec4D Zero = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D One = new MVec4D(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1); + + public static readonly MVec4D Unit = new MVec4D(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitX = new MVec4D(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitY = new MVec4D(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitZ = new MVec4D(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitW = new MVec4D(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitXy = new MVec4D(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitXz = new MVec4D(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitXw = new MVec4D(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitYz = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitYw = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitZw = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0); + public static readonly MVec4D UnitXyz = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0); + public static readonly MVec4D UnitXyw = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0); + public static readonly MVec4D UnitXzw = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0); + public static readonly MVec4D UnitYzw = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0); + public static readonly MVec4D UnitXyzw = new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1); + + #endregion + + public double S, X, Y, Z, W, Xy, Xz, Xw, Yz, Yw, Zw, Xyz, Xyw, Xzw, Yzw, Xyzw; + + public double Norm2 => S * S + X * X + Y * Y + Z * Z + W * W + Xy * Xy + Xz * Xz + Xw * Xw + Yz * Yz + Yw * Yw + Zw * Zw + Xyz * Xyz + + Xyw * Xyw + Xzw * Xzw + Yzw * Yzw + Xyzw * Xyzw; + + public double Norm => Math.Sqrt(Norm2); + + #region Swizzle + + public double Yx + { + get => -Xy; + set => Xy = -value; + } + + public double Zx + { + get => -Xz; + set => Xz = -value; + } + + public double Wx + { + get => -Xw; + set => Xw = -value; + } + + public double Zy + { + get => -Yz; + set => Yz = -value; + } + + public double Wy + { + get => -Yw; + set => Yw = -value; + } + + public double Wz + { + get => -Zw; + set => Zw = -value; + } + + public double Xzy + { + get => -Xyz; + set => Xyz = -value; + } + + public double Yxz + { + get => -Xyz; + set => Xyz = -value; + } + + public double Yzx + { + get => Xyz; + set => Xyz = value; + } + + public double Zxy + { + get => Xyz; + set => Xyz = value; + } + + public double Zyx + { + get => -Xyz; + set => Xyz = -value; + } + + public double Xwy + { + get => -Xyw; + set => Xyw = -value; + } + + public double Yxw + { + get => -Xyw; + set => Xyw = -value; + } + + public double Ywx + { + get => Xyw; + set => Xyw = value; + } + + public double Wxy + { + get => Xyw; + set => Xyw = value; + } + + public double Wyx + { + get => -Xyw; + set => Xyw = -value; + } + + public double Xwz + { + get => -Xzw; + set => Xzw = -value; + } + + public double Zxw + { + get => -Xzw; + set => Xzw = -value; + } + + public double Zwx + { + get => Xzw; + set => Xzw = value; + } + + public double Wxz + { + get => Xzw; + set => Xzw = value; + } + + public double Wzx + { + get => -Xzw; + set => Xzw = -value; + } + + public double Ywz + { + get => -Yzw; + set => Yzw = -value; + } + + public double Zyw + { + get => -Yzw; + set => Yzw = -value; + } + + public double Zwy + { + get => Yzw; + set => Yzw = value; + } + + public double Wyz + { + get => Yzw; + set => Yzw = value; + } + + public double Wzy + { + get => -Yzw; + set => Yzw = -value; + } + + public double Xywz + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Xzyw + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Xzwy + { + get => Xyzw; + set => Xyzw = value; + } + + public double Xwyz + { + get => Xyzw; + set => Xyzw = value; + } + + public double Xwzy + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Yxzw + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Yxwz + { + get => Xyzw; + set => Xyzw = value; + } + + public double Yzxw + { + get => Xyzw; + set => Xyzw = value; + } + + public double Yzwx + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Ywxz + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Ywzx + { + get => Xyzw; + set => Xyzw = value; + } + + public double Zxyw + { + get => Xyzw; + set => Xyzw = value; + } + + public double Zxwy + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Zyxw + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Zywx + { + get => Xyzw; + set => Xyzw = value; + } + + public double Zwxy + { + get => Xyzw; + set => Xyzw = value; + } + + public double Zwyx + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Wxyz + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Wxzy + { + get => Xyzw; + set => Xyzw = value; + } + + public double Wyxz + { + get => Xyzw; + set => Xyzw = value; + } + + public double Wyzx + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Wzxy + { + get => -Xyzw; + set => Xyzw = -value; + } + + public double Wzyx + { + get => Xyzw; + set => Xyzw = value; + } + + #endregion + + #region Inverses + + public MVec4D Neg => new MVec4D(-S, -X, -Y, -Z, -W, -Xy, -Xz, -Xw, -Yz, -Yw, -Zw, -Xyz, -Xyw, -Xzw, -Yzw, -Xyzw); + + public MVec4D? Inv + { + get + { + var v = Mul(this, Rev); + if (v == Vec0(v.S)) + return Mul(1d / v.S, Rev); + + return null; + } + } + + public MVec4D Rev => new MVec4D(S, X, Y, Z, W, Yx, Zx, Wx, Zy, Wy, Wz, Zyx, Wyx, Wzx, Wzy, Wzyx); + + public MVec4D Conj => new MVec4D(Xyzw, -Yzw, Xzw, -Xyw, -Xyz, Zw, Yw, Yz, Xw, Xz, Xy, -W, Z, Y, X, S); + + public MVec4D Dual => new MVec4D(Xyzw, -Yzw, -Xzw, -Xyw, Xyz, -Zw, Yw, -Yz, -Xw, Xz, -Xy, W, -Z, Y, -Xz, S); + + #endregion + + #region Constructor + + public MVec4D(double s, double x, double y, double z, double w, double xy, double xz, double xw, double yz, double yw, double zw, double xyz, + double xyw, double xzw, double yzw, double xyzw) + { + S = s; + X = x; + Y = y; + Z = z; + W = w; + Xy = xy; + Xz = xz; + Xw = xw; + Yz = yz; + Yw = yw; + Zw = zw; + Xyz = xyz; + Xyw = xyw; + Xzw = xzw; + Yzw = yzw; + Xyzw = xyzw; + } + + public MVec4D(MVec4D s, MVec4D v, MVec4D b, MVec4D t, MVec4D q) + { + S = s.S; + X = v.X; + Y = v.Y; + Z = v.Z; + W = v.W; + Xy = b.Xy; + Xz = b.Xz; + Xw = b.Xw; + Yz = b.Yz; + Yw = b.Yw; + Zw = b.Zw; + Xyz = t.Xyz; + Xyw = t.Xyw; + Xzw = t.Xzw; + Yzw = t.Yzw; + Xyzw = q.Xyzw; + } + + #endregion + + #region Factory Methods + + public static MVec4D Vec0(double s) + { + return new MVec4D(s, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + } + + public static MVec4D Vec1(double x, double y, double z, double w) + { + return new MVec4D(0, x, y, z, w, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); + } + + public static MVec4D Vec2(double xy, double xz, double xw, double yz, double yw, double zw) + { + return new MVec4D(0, 0, 0, 0, 0, xy, xz, xw, yz, yw, zw, 0, 0, 0, 0, 0); + } + + public static MVec4D Vec3(double xyz, double xyw, double xzw, double yzw) + { + return new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, xyz, xyw, xzw, yzw, 0); + } + + public static MVec4D Vec4(double xyzw) + { + return new MVec4D(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, xyzw); + } + + public static MVec4D Rotor(double angle, MVec4D plane) + { + return Add(Math.Cos(angle), Mul(Math.Sin(angle), plane)); + } + + #endregion + + #region Operations + + public static MVec4D Add(MVec4D m, MVec4D n) + { + return new MVec4D( + m.S + n.S, + m.X + n.X, + m.Y + n.Y, + m.Z + n.Z, + m.W + n.W, + m.Xy + n.Xy, + m.Xz + n.Xz, + m.Xw + n.Xw, + m.Yz + n.Yz, + m.Yw + n.Yw, + m.Zw + n.Zw, + m.Xyz + n.Xyz, + m.Xyw + n.Xyw, + m.Xzw + n.Xzw, + m.Yzw + n.Yzw, + m.Xyzw + n.Xyzw + ); + } + + public static MVec4D Add(MVec4D m, double c) => Add(c, m); + + public static MVec4D Add(double c, MVec4D m) + { + return new MVec4D( + m.S + c, + m.X, + m.Y, + m.Z, + m.W, + m.Xy, + m.Xz, + m.Xw, + m.Yz, + m.Yw, + m.Zw, + m.Xyz, + m.Xyw, + m.Xzw, + m.Yzw, + m.Xyzw + ); + } + + public static MVec4D Sub(MVec4D m, MVec4D n) + { + return new MVec4D( + m.S - n.S, + m.X - n.X, + m.Y - n.Y, + m.Z - n.Z, + m.W - n.W, + m.Xy - n.Xy, + m.Xz - n.Xz, + m.Xw - n.Xw, + m.Yz - n.Yz, + m.Yw - n.Yw, + m.Zw - n.Zw, + m.Xyz - n.Xyz, + m.Xyw - n.Xyw, + m.Xzw - n.Xzw, + m.Yzw - n.Yzw, + m.Xyzw - n.Xyzw + ); + } + + public static MVec4D Sub(MVec4D m, double c) + { + return new MVec4D( + m.S - c, + m.X, + m.Y, + m.Z, + m.W, + m.Xy, + m.Xz, + m.Xw, + m.Yz, + m.Yw, + m.Zw, + m.Xyz, + m.Xyw, + m.Xzw, + m.Yzw, + m.Xyzw + ); + } + + public static MVec4D Sub(double c, MVec4D m) + { + return new MVec4D( + c - m.S, + -m.X, + -m.Y, + -m.Z, + -m.W, + -m.Xy, + -m.Xz, + -m.Xw, + -m.Yz, + -m.Yw, + -m.Zw, + -m.Xyz, + -m.Xyw, + -m.Xzw, + -m.Yzw, + -m.Xyzw + ); + } + + public static MVec4D Mul(MVec4D m, MVec4D n) + { + var r = new MVec4D(); + if (Math.Abs(m.S) > 1e-10) + { + r.S += m.S * n.S; + r.X += m.S * n.X; + r.Y += m.S * n.Y; + r.Z += m.S * n.Z; + r.W += m.S * n.W; + r.Xy += m.S * n.Xy; + r.Xz += m.S * n.Xz; + r.Xw += m.S * n.Xw; + r.Yz += m.S * n.Yz; + r.Yw += m.S * n.Yw; + r.Zw += m.S * n.Zw; + r.Xyz += m.S * n.Xyz; + r.Xyw += m.S * n.Xyw; + r.Xzw += m.S * n.Xzw; + r.Yzw += m.S * n.Yzw; + r.Xyzw += m.S * n.Xyzw; + } + + if (Math.Abs(m.X) > 1e-10) + { + r.S += m.X * n.X; + r.X += m.X * n.S; + r.Y += m.X * n.Xy; + r.Z += m.X * n.Xz; + r.W += m.X * n.Xw; + r.Xy += m.X * n.Y; + r.Xz += m.X * n.Z; + r.Xw += m.X * n.W; + r.Yz += m.X * n.Xyz; + r.Yw += m.X * n.Xyw; + r.Zw += m.X * n.Xzw; + r.Xyz += m.X * n.Yz; + r.Xyw += m.X * n.Yw; + r.Xzw += m.X * n.Zw; + r.Yzw += m.X * n.Xyzw; + r.Xyzw += m.X * n.Yzw; + } + + if (Math.Abs(m.Y) > 1e-10) + { + r.S += m.Y * n.Y; + r.X += m.Y * n.Yx; + r.Y += m.Y * n.S; + r.Z += m.Y * n.Yz; + r.W += m.Y * n.Yw; + r.Yx += m.Y * n.X; + r.Xz += m.Y * n.Yxz; + r.Xw += m.Y * n.Yxw; + r.Yz += m.Y * n.Z; + r.Yw += m.Y * n.W; + r.Zw += m.Y * n.Yzw; + r.Yxz += m.Y * n.Xz; + r.Yxw += m.Y * n.Xw; + r.Xzw += m.Y * n.Yxzw; + r.Yzw += m.Y * n.Zw; + r.Yxzw += m.Y * n.Xzw; + } + + if (Math.Abs(m.Z) > 1e-10) + { + r.S += m.Z * n.Z; + r.X += m.Z * n.Zx; + r.Y += m.Z * n.Zy; + r.Z += m.Z * n.S; + r.W += m.Z * n.Zw; + r.Xy += m.Z * n.Zxy; + r.Zx += m.Z * n.X; + r.Xw += m.Z * n.Zxw; + r.Zy += m.Z * n.Y; + r.Yw += m.Z * n.Zyw; + r.Zw += m.Z * n.W; + r.Zxy += m.Z * n.Xy; + r.Xyw += m.Z * n.Zxyw; + r.Zxw += m.Z * n.Xw; + r.Zyw += m.Z * n.Yw; + r.Zxyw += m.Z * n.Xyw; + } + + if (Math.Abs(m.W) > 1e-10) + { + r.S += m.W * n.W; + r.X += m.W * n.Wx; + r.Y += m.W * n.Wy; + r.Z += m.W * n.Wz; + r.W += m.W * n.S; + r.Xy += m.W * n.Wxy; + r.Xz += m.W * n.Wxz; + r.Wx += m.W * n.X; + r.Yz += m.W * n.Wyz; + r.Wy += m.W * n.Y; + r.Wz += m.W * n.Z; + r.Xyz += m.W * n.Wxyz; + r.Wxy += m.W * n.Xy; + r.Wxz += m.W * n.Xz; + r.Wyz += m.W * n.Yz; + r.Wxyz += m.W * n.Xyz; + } + + if (Math.Abs(m.Xy) > 1e-10) + { + r.S += m.Xy * n.Yx; + r.X += m.Xy * n.Y; + r.Y += m.Yx * n.X; + r.Z += m.Xy * n.Yxz; + r.W += m.Xy * n.Yxw; + r.Xy += m.Xy * n.S; + r.Xz += m.Xy * n.Yz; + r.Xw += m.Xy * n.Yw; + r.Yz += m.Yx * n.Xz; + r.Yw += m.Yx * n.Xw; + r.Zw += m.Xy * n.Yxzw; + r.Xyz += m.Xy * n.Z; + r.Xyw += m.Xy * n.W; + r.Xzw += m.Xy * n.Yzw; + r.Yzw += m.Yx * n.Xzw; + r.Xyzw += m.Xy * n.Zw; + } + + if (Math.Abs(m.Xz) > 1e-10) + { + r.S += m.Xz * n.Zx; + r.X += m.Xz * n.Z; + r.Y += m.Xz * n.Zxy; + r.Z += m.Zx * n.X; + r.W += m.Xz * n.Zxw; + r.Xy += m.Xz * n.Zy; + r.Xz += m.Xz * n.S; + r.Xw += m.Xz * n.Zw; + r.Zy += m.Zx * n.Xy; + r.Yw += m.Xz * n.Zxyw; + r.Zw += m.Zx * n.Xw; + r.Xzy += m.Xz * n.Y; + r.Xyw += m.Xz * n.Zyw; + r.Xzw += m.Xz * n.W; + r.Zyw += m.Zx * n.Xyw; + r.Xzyw += m.Xz * n.Yw; + } + + if (Math.Abs(m.Xw) > 1e-10) + { + r.S += m.Xw * n.Wx; + r.X += m.Xw * n.W; + r.Y += m.Xw * n.Wxy; + r.Z += m.Xw * n.Wxz; + r.W += m.Wx * n.X; + r.Xy += m.Xw * n.Wy; + r.Xz += m.Xw * n.Wz; + r.Xw += m.Xw * n.S; + r.Yz += m.Xw * n.Wxyz; + r.Wy += m.Wx * n.Xy; + r.Wz += m.Wx * n.Xz; + r.Xyz += m.Xw * n.Wyz; + r.Xwy += m.Xw * n.Y; + r.Xwz += m.Xw * n.Z; + r.Wyz += m.Wx * n.Xyz; + r.Xwyz += m.Xw * n.Yz; + } + + if (Math.Abs(m.Yz) > 1e-10) + { + r.S += m.Yz * n.Zy; + r.X += m.Yz * n.Zyx; + r.Y += m.Yz * n.Z; + r.Z += m.Zy * n.Y; + r.W += m.Yz * n.Zyw; + r.Yx += m.Yz * n.Zx; + r.Zx += m.Zy * n.Yx; + r.Xw += m.Yz * n.Zyxw; + r.Yz += m.Yz * n.S; + r.Yw += m.Yz * n.Zw; + r.Zw += m.Zy * n.Yw; + r.Yzx += m.Yz * n.X; + r.Yxw += m.Yz * n.Zxw; + r.Zxw += m.Zy * n.Yxw; + r.Yzw += m.Yz * n.W; + r.Yzxw += m.Yz * n.Xw; + } + + if (Math.Abs(m.Yw) > 1e-10) + { + r.S += m.Yw * n.Wy; + r.X += m.Yw * n.Wyx; + r.Y += m.Yw * n.W; + r.Z += m.Yw * n.Wyz; + r.W += m.Wy * n.Y; + r.Yx += m.Yw * n.Wx; + r.Xz += m.Yw * n.Wyxz; + r.Wx += m.Wy * n.Yx; + r.Yz += m.Yw * n.Wz; + r.Yw += m.Yw * n.S; + r.Wz += m.Wy * n.Yz; + r.Yxz += m.Yw * n.Wxz; + r.Ywx += m.Yw * n.X; + r.Wxz += m.Wy * n.Yxz; + r.Ywz += m.Yw * n.Z; + r.Ywxz += m.Yw * n.Xz; + } + + if (Math.Abs(m.Zw) > 1e-10) + { + r.S += m.Zw * n.Wz; + r.X += m.Zw * n.Wzx; + r.Y += m.Zw * n.Wzy; + r.Z += m.Zw * n.W; + r.W += m.Wz * n.Z; + r.Xy += m.Zw * n.Wzxy; + r.Zx += m.Zw * n.Wx; + r.Wx += m.Wz * n.Zx; + r.Zy += m.Zw * n.Wy; + r.Wy += m.Wz * n.Zy; + r.Zw += m.Zw * n.S; + r.Zxy += m.Zw * n.Wxy; + r.Wxy += m.Wz * n.Zxy; + r.Zwx += m.Zw * n.X; + r.Zwy += m.Zw * n.Y; + r.Zwxy += m.Zw * n.Xy; + } + + if (Math.Abs(m.Xyz) > 1e-10) + { + r.S += m.Xyz * n.Zyx; + r.X += m.Xyz * n.Zy; + r.Y += m.Yxz * n.Zx; + r.Z += m.Zxy * n.Yx; + r.W += m.Xyz * n.Zyxw; + r.Xy += m.Xyz * n.Z; + r.Xz += m.Xzy * n.Y; + r.Xw += m.Xyz * n.Zyw; + r.Yz += m.Yzx * n.X; + r.Yw += m.Yxz * n.Zxw; + r.Zw += m.Zxy * n.Yxw; + r.Xyz += m.Xyz * n.S; + r.Xyw += m.Xyz * n.Zw; + r.Xzw += m.Xzy * n.Yw; + r.Yzw += m.Yzx * n.Xw; + r.Xyzw += m.Xyz * n.W; + } + + if (Math.Abs(m.Xyw) > 1e-10) + { + r.S += m.Xyw * n.Wyx; + r.X += m.Xyw * n.Wy; + r.Y += m.Yxw * n.Wx; + r.Z += m.Xyw * n.Wyxz; + r.W += m.Wxy * n.Yx; + r.Xy += m.Xyw * n.W; + r.Xz += m.Xyw * n.Wyz; + r.Xw += m.Xwy * n.Y; + r.Yz += m.Yxw * n.Wxz; + r.Yw += m.Ywx * n.X; + r.Wz += m.Wxy * n.Yxz; + r.Xyz += m.Xyw * n.Wz; + r.Xyw += m.Xyw * n.S; + r.Xwz += m.Xwy * n.Yz; + r.Ywz += m.Ywx * n.Xz; + r.Xywz += m.Xyw * n.Z; + } + + if (Math.Abs(m.Xzw) > 1e-10) + { + r.S += m.Xzw * n.Wzx; + r.X += m.Xzw * n.Wz; + r.Y += m.Xzw * n.Wzxy; + r.Z += m.Zxw * n.Wx; + r.W += m.Wxz * n.Zx; + r.Xy += m.Xzw * n.Wzy; + r.Xz += m.Xzw * n.W; + r.Xw += m.Xwz * n.Z; + r.Zy += m.Zxw * n.Wxy; + r.Wy += m.Wxz * n.Zxy; + r.Zw += m.Zwx * n.X; + r.Xzy += m.Xzw * n.Wy; + r.Xwy += m.Xwz * n.Zy; + r.Xzw += m.Xzw * n.S; + r.Zwy += m.Zwx * n.Xy; + r.Xzwy += m.Xzw * n.Y; + } + + if (Math.Abs(m.Yzw) > 1e-10) + { + r.S += m.Yzw * n.Wzy; + r.X += m.Yzw * n.Wzyx; + r.Y += m.Yzw * n.Wz; + r.Z += m.Zyw * n.Wy; + r.W += m.Wyz * n.Zy; + r.Yx += m.Yzw * n.Wzx; + r.Zx += m.Zyw * n.Wyx; + r.Wx += m.Wyz * n.Zyx; + r.Yz += m.Yzw * n.W; + r.Yw += m.Ywz * n.Z; + r.Zw += m.Zwy * n.Y; + r.Yzx += m.Yzw * n.Wx; + r.Ywx += m.Ywz * n.Zx; + r.Zwx += m.Zwy * n.Yx; + r.Yzw += m.Yzw * n.S; + r.Yzwx += m.Yzw * n.X; + } + + if (Math.Abs(m.Xyzw) > 1e-10) + { + r.S += m.Xyzw * n.Wzyx; + r.X += m.Xyzw * n.Wzy; + r.Y += m.Yxzw * n.Wzx; + r.Z += m.Zxyw * n.Wyx; + r.W += m.Wxyz * n.Zyx; + r.Xy += m.Xyzw * n.Wz; + r.Xz += m.Xzyw * n.Wy; + r.Xw += m.Xwyz * n.Zy; + r.Yz += m.Yzxw * n.Wx; + r.Yw += m.Ywxz * n.Zx; + r.Zw += m.Zwxy * n.Yx; + r.Xyz += m.Xyzw * n.W; + r.Xyw += m.Xywz * n.Z; + r.Xzw += m.Xzwy * n.Y; + r.Yzw += m.Yzwx * n.X; + r.Xyzw += m.Xyzw * n.S; + } + + return r; + } + + public static MVec4D Mul(MVec4D m, double c) => Mul(c, m); + + public static MVec4D Mul(double c, MVec4D m) + { + return new MVec4D( + c * m.S, + c * m.X, + c * m.Y, + c * m.Z, + c * m.W, + c * m.Xy, + c * m.Xz, + c * m.Xw, + c * m.Yz, + c * m.Yw, + c * m.Zw, + c * m.Xyz, + c * m.Xyw, + c * m.Xzw, + c * m.Yzw, + c * m.Xyzw + ); + } + + public static MVec4D? Div(MVec4D m, MVec4D n) + { + var inv = n.Inv; + if (!inv.HasValue) + return null; + + return Mul(m, inv.Value); + } + + public static MVec4D Div(MVec4D m, double c) => Mul(m, 1 / c); + + public static MVec4D? Div(double c, MVec4D m) + { + var inv = m.Inv; + if (!inv.HasValue) + return null; + return Mul(c, inv.Value); + } + + #endregion + + #region Equality + + public static bool operator ==(MVec4D m, MVec4D n) + { + var diff = Sub(m, n).Norm2; + return Math.Abs(diff) <= 10e-8; + } + + public static bool operator !=(MVec4D m, MVec4D n) + { + var diff = Sub(m, n).Norm2; + return Math.Abs(diff) > 10e-8; + } + + #endregion + + #region Formatting + + public override string ToString() + { + var sb = new List(); + if (Math.Abs(S) > 10e-2) sb.Add($"{S:F2}"); + if (Math.Abs(X) > 10e-2) sb.Add($"{X:F2}e1"); + if (Math.Abs(Y) > 10e-2) sb.Add($"{Y:F2}e2"); + if (Math.Abs(Z) > 10e-2) sb.Add($"{Z:F2}e3"); + if (Math.Abs(W) > 10e-2) sb.Add($"{W:F2}e4"); + if (Math.Abs(Xy) > 10e-2) sb.Add($"{Xy:F2}e12"); + if (Math.Abs(Yz) > 10e-2) sb.Add($"{Yz:F2}e23"); + if (Math.Abs(Zw) > 10e-2) sb.Add($"{Zw:F2}e34"); + if (Math.Abs(Wx) > 10e-2) sb.Add($"{Wx:F2}e41"); + if (Math.Abs(Xz) > 10e-2) sb.Add($"{Xz:F2}e13"); + if (Math.Abs(Yw) > 10e-2) sb.Add($"{Yw:F2}e24"); + if (Math.Abs(Xyz) > 10e-2) sb.Add($"{Xyz:F2}e123"); + if (Math.Abs(Yzw) > 10e-2) sb.Add($"{Yzw:F2}e234"); + if (Math.Abs(Zwx) > 10e-2) sb.Add($"{Zwx:F2}e341"); + if (Math.Abs(Wxy) > 10e-2) sb.Add($"{Wxy:F2}e412"); + if (Math.Abs(Xyzw) > 10e-2) sb.Add($"{Xyzw:F2}e1234"); + return $"({string.Join(" + ", sb)})"; + } + + #endregion + } +} \ No newline at end of file diff --git a/Tetrahedrons/Program.cs b/Tetrahedrons/Program.cs index 6175764..2015418 100644 --- a/Tetrahedrons/Program.cs +++ b/Tetrahedrons/Program.cs @@ -1,13 +1,28 @@ -using OpenTK; +using System; +using OpenTK; namespace Tetrahedrons { - internal class Program : GameWindow - { - public static void Main(string[] args) - { - using(var p = new TetrahedronWindow()) - p.Run(); - } - } + internal class Program : GameWindow + { + public static void NotMain(string[] args) + { + using (var p = new TetrahedronWindow()) + p.Run(); + } + + public static void Main(string[] args) + { + var m = new MVec4D(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0); + var i = m.Inv ?? MVec4D.Zero; + + Console.Out.WriteLine("m = {0}", m); + Console.Out.WriteLine("m.Rev = {0}", m.Rev); + Console.Out.WriteLine("MVec4D.Mul(m, m.Rev) = {0}", MVec4D.Mul(m, m.Rev)); + Console.Out.WriteLine("m.Inv = {0}", i); + + Console.Out.WriteLine("MVec4D.Mul(m, m.Inv) = {0}", MVec4D.Mul(m, i)); + Console.Out.WriteLine("MVec4D.Mul(m.Inv, m) = {0}", MVec4D.Mul(i, m)); + } + } } \ No newline at end of file diff --git a/Tetrahedrons/TetrahedronWindow.cs b/Tetrahedrons/TetrahedronWindow.cs index 1ea9e16..d7061ed 100644 --- a/Tetrahedrons/TetrahedronWindow.cs +++ b/Tetrahedrons/TetrahedronWindow.cs @@ -34,7 +34,7 @@ namespace Tetrahedrons GL.End(); } - public static Polygon operator +(Polygon p, MVec3d m) + public static Polygon operator +(Polygon p, MVec3D m) { var pts = new Vector3d[p.Points.Length]; for (var i = 0; i < p.Points.Length; i++) @@ -42,7 +42,7 @@ namespace Tetrahedrons return new Polygon(pts, p.Color); } - public static Polygon operator -(Polygon p, MVec3d m) + public static Polygon operator -(Polygon p, MVec3D m) { var pts = new Vector3d[p.Points.Length]; for (var i = 0; i < p.Points.Length; i++) @@ -50,7 +50,7 @@ namespace Tetrahedrons return new Polygon(pts, p.Color); } - public static Polygon operator *(Polygon p, MVec3d m) + public static Polygon operator *(Polygon p, MVec3D m) { var pts = new Vector3d[p.Points.Length]; for (var i = 0; i < p.Points.Length; i++) @@ -58,7 +58,7 @@ namespace Tetrahedrons return new Polygon(pts, p.Color); } - public static Polygon operator ^(Polygon p, MVec3d m) + public static Polygon operator ^(Polygon p, MVec3D m) { var pts = new Vector3d[p.Points.Length]; for (var i = 0; i < p.Points.Length; i++) @@ -66,7 +66,7 @@ namespace Tetrahedrons return new Polygon(pts, p.Color); } - public static Polygon operator &(Polygon p, MVec3d m) + public static Polygon operator &(Polygon p, MVec3D m) { var pts = new Vector3d[p.Points.Length]; for (var i = 0; i < p.Points.Length; i++) @@ -130,7 +130,7 @@ namespace Tetrahedrons } } - public Polygon Intersection(MVec3d blade, Vector3d pivot) + public Polygon Intersection(MVec3D blade, Vector3d pivot) { var pts = new Vector3d[4]; for (var i = 0; i < 4; i++) @@ -138,7 +138,7 @@ namespace Tetrahedrons pts[i] = Points[i] - pivot; } - var a = blade * MVec3d.Unit123; + var a = blade * MVec3D.Unit123; var sides = new double[4]; for (var i = 0; i < 4; i++) { @@ -175,9 +175,9 @@ namespace Tetrahedrons private Tetrahedron _tetra; private Polygon _polyg; - private MVec3d _b1 = MVec3d.Unit1; - private MVec3d _b2 = MVec3d.Unit3; - private MVec3d _pivot = MVec3d.Zero; + private MVec3D _b1 = MVec3D.Unit1; + private MVec3D _b2 = MVec3D.Unit3; + private MVec3D _pivot = MVec3D.Zero; private Polygon _square = new Polygon(new[] {new Vector3d(-1, -1, 0), new Vector3d(-1, 1, 0), new Vector3d(1, -1, 0), new Vector3d(1, 1, 0),}, Color.Red); @@ -232,7 +232,7 @@ namespace Tetrahedrons GL.Begin(PrimitiveType.Lines); GL.Vertex3(_pivot.VectorPart); - GL.Vertex3((_pivot + (_b2 ^ _b1) * MVec3d.Unit123 * .5).VectorPart); + GL.Vertex3((_pivot + (_b2 ^ _b1) * MVec3D.Unit123 * .5).VectorPart); GL.End(); GL.Enable(EnableCap.DepthTest); @@ -255,26 +255,26 @@ namespace Tetrahedrons if (Keyboard[Key.Left]) { - var rot = MVec3d.Complex(-e.Time, MVec3d.Unit12); + var rot = MVec3D.Complex(-e.Time, MVec3D.Unit12); _b1 = rot * _b1 * ~rot; // todo make blade classes _b2 = rot * _b2 * ~rot; } if (Keyboard[Key.Right]) { - var rot = MVec3d.Complex(e.Time, MVec3d.Unit12); + var rot = MVec3D.Complex(e.Time, MVec3D.Unit12); _b1 = rot * _b1 * ~rot; _b2 = rot * _b2 * ~rot; } if (Keyboard[Key.Up]) { - _pivot -= (_b1 ^ _b2) * MVec3d.Unit123 * e.Time; + _pivot -= (_b1 ^ _b2) * MVec3D.Unit123 * e.Time; } if (Keyboard[Key.Down]) { - _pivot += (_b1 ^ _b2) * MVec3d.Unit123 * e.Time; + _pivot += (_b1 ^ _b2) * MVec3D.Unit123 * e.Time; } if (_pause) @@ -310,7 +310,7 @@ namespace Tetrahedrons for (var i = 0; i < _polyg.Points.Length; i++) { - MVec3d p = _polyg.Points[i]; + MVec3D p = _polyg.Points[i]; var x = ((p & _b1) * ~_b1) & _b1; var y = ((p ^ _b1) * ~_b1) & ((_b2 ^ _b1) * ~_b1); _polyg.Points[i] = new Vector3d(x.E - px.E, y.E - py.E, 0); @@ -331,9 +331,9 @@ namespace Tetrahedrons _div = !_div; break; case Key.R: - _b1 = MVec3d.Unit1; - _b2 = MVec3d.Unit3; - _pivot = MVec3d.Zero; + _b1 = MVec3D.Unit1; + _b2 = MVec3D.Unit3; + _pivot = MVec3D.Zero; break; } } diff --git a/Tetrahedrons/Tetrahedrons.csproj b/Tetrahedrons/Tetrahedrons.csproj index a76912a..70b1a79 100644 --- a/Tetrahedrons/Tetrahedrons.csproj +++ b/Tetrahedrons/Tetrahedrons.csproj @@ -43,6 +43,7 @@ + From 8bec0d6b805669510f03133a87230d269882f0a4 Mon Sep 17 00:00:00 2001 From: allem Date: Tue, 3 Apr 2018 23:40:06 -0400 Subject: [PATCH 2/7] inner and outer products --- Tetrahedrons/MVec3d.cs | 12 +- Tetrahedrons/MVec4D.cs | 177 +++++++++++++++++++++++++++--- Tetrahedrons/TetrahedronWindow.cs | 4 +- 3 files changed, 166 insertions(+), 27 deletions(-) diff --git a/Tetrahedrons/MVec3d.cs b/Tetrahedrons/MVec3d.cs index 4ce9eb8..e4c2cf9 100644 --- a/Tetrahedrons/MVec3d.cs +++ b/Tetrahedrons/MVec3d.cs @@ -68,17 +68,7 @@ namespace Tetrahedrons return new MVec3D(0, 0, 0, 0, 0, 0, 0, e123); } - public static MVec3D Complex(double real, double imag) - { - return new MVec3D(real, 0, 0, 0, 0, 0, 0, imag); - } - - public static MVec3D Complex(double w, double i, double j, double k) - { - return new MVec3D(w, 0, 0, 0, i, j, k, 0); - } - - public static MVec3D Complex(double angle, MVec3D plane) + public static MVec3D Rotor(double angle, MVec3D plane) { return Math.Cos(angle) + Math.Sin(angle) * plane; } diff --git a/Tetrahedrons/MVec4D.cs b/Tetrahedrons/MVec4D.cs index 62216a0..ddc59be 100644 --- a/Tetrahedrons/MVec4D.cs +++ b/Tetrahedrons/MVec4D.cs @@ -1,8 +1,11 @@ using System; using System.Collections.Generic; +using System.Data.SqlClient; +using System.Dynamic; using System.Runtime.Remoting.Messaging; using System.Text; using OpenTK; +using OpenTK.Graphics.ES30; namespace Tetrahedrons { @@ -38,7 +41,26 @@ namespace Tetrahedrons Xyw * Xyw + Xzw * Xzw + Yzw * Yzw + Xyzw * Xyzw; public double Norm => Math.Sqrt(Norm2); - + + public MVec4D Grade(int k) + { + switch (k) + { + case 0: + return Vec0(S); + case 1: + return Vec1(X, Y, Z, W); + case 2: + return Vec2(Xy, Xz, Xw, Yz, Yw, Zw); + case 3: + return Vec3(Xyz, Xyw, Xzw, Yzw); + case 4: + return Vec4(Xyzw); + default: + return Zero; + } + } + #region Swizzle public double Yx @@ -904,25 +926,152 @@ namespace Tetrahedrons ); } - public static MVec4D? Div(MVec4D m, MVec4D n) - { - var inv = n.Inv; - if (!inv.HasValue) - return null; - - return Mul(m, inv.Value); - } + public static MVec4D? Div(MVec4D m, MVec4D n) => Mul(m, n.Inv); public static MVec4D Div(MVec4D m, double c) => Mul(m, 1 / c); - public static MVec4D? Div(double c, MVec4D m) + public static MVec4D? Div(double c, MVec4D n) => Mul(c, n.Inv); + + public static MVec4D Inner(MVec4D m, MVec4D n) { - var inv = m.Inv; - if (!inv.HasValue) - return null; - return Mul(c, inv.Value); + var m0 = m.Grade(0); + var m1 = m.Grade(1); + var m2 = m.Grade(2); + var m3 = m.Grade(3); + var m4 = m.Grade(4); + + var n0 = n.Grade(0); + var n1 = n.Grade(1); + var n2 = n.Grade(2); + var n3 = n.Grade(3); + var n4 = n.Grade(4); + + return new MVec4D( + m0 * n0 + m1 * n1 + m2 * n2 + m3 * n3 + n4 * n4, + m0 * n1 + m1 * n2 + m2 * n3 + m3 * n4, + m0 * n2 + m1 * n3 + m2 * n4, + m0 * n3 + m1 * n4, + m0 * n4 + ); } + public static MVec4D Inner(MVec4D m, double c) => Mul(m, c); + + public static MVec4D Inner(double c, MVec4D n) => Mul(c, n); + + public static MVec4D Outer(MVec4D m, MVec4D n) + { + var m0 = m.Grade(0); + var m1 = m.Grade(1); + var m2 = m.Grade(2); + var m3 = m.Grade(3); + var m4 = m.Grade(4); + + var n0 = n.Grade(0); + var n1 = n.Grade(1); + var n2 = n.Grade(2); + var n3 = n.Grade(3); + var n4 = n.Grade(4); + + return new MVec4D( + m0 * n0, + m1 * n0 + m0 * n1, + m2 * n0 + m1 * n1 + m0 * n2, + m3 * n0 + m2 * n1 + m1 * n2 + m0 * n3, + m4 * n0 + m3 * n1 + m2 * n2 + m1 * n3 + m0 * n4 + ); + } + + public static MVec4D Outer(MVec4D m, double c) => Mul(m, c); + + public static MVec4D Outer(double c, MVec4D n) => Mul(c, n); + + + public static MVec4D? Add(MVec4D? m, MVec4D? n) => m.HasValue && n.HasValue ? (MVec4D?) Add(m.Value, n.Value) : null; + public static MVec4D? Add(MVec4D? m, double? n) => m.HasValue && n.HasValue ? (MVec4D?) Add(m.Value, n.Value) : null; + public static MVec4D? Add(double? m, MVec4D? n) => n.HasValue && m.HasValue ? (MVec4D?) Add(m.Value, n.Value) : null; + + public static MVec4D? Sub(MVec4D? m, MVec4D? n) => m.HasValue && n.HasValue ? (MVec4D?) Sub(m.Value, n.Value) : null; + public static MVec4D? Sub(MVec4D? m, double? n) => m.HasValue && n.HasValue ? (MVec4D?) Sub(m.Value, n.Value) : null; + public static MVec4D? Sub(double? m, MVec4D? n) => n.HasValue && m.HasValue ? (MVec4D?) Sub(m.Value, n.Value) : null; + + public static MVec4D? Mul(MVec4D? m, MVec4D? n) => m.HasValue && n.HasValue ? (MVec4D?) Mul(m.Value, n.Value) : null; + public static MVec4D? Mul(MVec4D? m, double? n) => m.HasValue && n.HasValue ? (MVec4D?) Mul(m.Value, n.Value) : null; + public static MVec4D? Mul(double? m, MVec4D? n) => n.HasValue && m.HasValue ? (MVec4D?) Mul(m.Value, n.Value) : null; + + public static MVec4D? Div(MVec4D? m, MVec4D? n) => m.HasValue && n.HasValue ? (MVec4D?) Div(m.Value, n.Value) : null; + public static MVec4D? Div(MVec4D? m, double? n) => m.HasValue && n.HasValue ? (MVec4D?) Div(m.Value, n.Value) : null; + public static MVec4D? Div(double? m, MVec4D? n) => n.HasValue && m.HasValue ? (MVec4D?) Div(m.Value, n.Value) : null; + + public static MVec4D? Inner(MVec4D? m, MVec4D? n) => m.HasValue && n.HasValue ? (MVec4D?) Inner(m.Value, n.Value) : null; + public static MVec4D? Inner(MVec4D? m, double? n) => m.HasValue && n.HasValue ? (MVec4D?) Inner(m.Value, n.Value) : null; + public static MVec4D? Inner(double? m, MVec4D? n) => n.HasValue && m.HasValue ? (MVec4D?) Inner(m.Value, n.Value) : null; + + public static MVec4D? Outer(MVec4D? m, MVec4D? n) => m.HasValue && n.HasValue ? (MVec4D?) Outer(m.Value, n.Value) : null; + public static MVec4D? Outer(MVec4D? m, double? n) => m.HasValue && n.HasValue ? (MVec4D?) Outer(m.Value, n.Value) : null; + public static MVec4D? Outer(double? m, MVec4D? n) => n.HasValue && m.HasValue ? (MVec4D?) Outer(m.Value, n.Value) : null; + + #endregion + + #region Operators + + public static MVec4D operator +(MVec4D m, MVec4D n) => Add(m, n); + public static MVec4D operator +(double c, MVec4D n) => Add(c, n); + public static MVec4D operator +(MVec4D m, double c) => Add(m, c); + + public static MVec4D operator -(MVec4D m, MVec4D n) => Sub(m, n); + public static MVec4D operator -(double c, MVec4D n) => Sub(c, n); + public static MVec4D operator -(MVec4D m, double c) => Sub(m, c); + + public static MVec4D operator *(MVec4D m, MVec4D n) => Mul(m, n); + public static MVec4D operator *(double c, MVec4D n) => Mul(c, n); + public static MVec4D operator *(MVec4D m, double c) => Mul(m, c); + + public static MVec4D? operator /(MVec4D m, MVec4D n) => Div(m, n); + public static MVec4D? operator /(double c, MVec4D n) => Div(c, n); + public static MVec4D operator /(MVec4D m, double c) => Div(m, c); + + public static MVec4D operator &(MVec4D m, MVec4D n) => Inner(m, n); + public static MVec4D operator &(MVec4D m, double c) => Inner(m, c); + public static MVec4D operator &(double c, MVec4D n) => Inner(c, n); + + public static MVec4D operator ^(MVec4D m, MVec4D n) => Outer(m, n); + public static MVec4D operator ^(MVec4D m, double c) => Outer(m, c); + public static MVec4D operator ^(double c, MVec4D n) => Outer(c, n); + + public static MVec4D? operator +(MVec4D? m, MVec4D? n) => Add(m, n); + public static MVec4D? operator +(double? c, MVec4D? n) => Add(c, n); + public static MVec4D? operator +(MVec4D? m, double? c) => Add(m, c); + + public static MVec4D? operator -(MVec4D? m, MVec4D? n) => Sub(m, n); + public static MVec4D? operator -(double? c, MVec4D? n) => Sub(c, n); + public static MVec4D? operator -(MVec4D? m, double? c) => Sub(m, c); + + public static MVec4D? operator *(MVec4D? m, MVec4D? n) => Mul(m, n); + public static MVec4D? operator *(double? c, MVec4D? n) => Mul(c, n); + public static MVec4D? operator *(MVec4D? m, double? c) => Mul(m, c); + + public static MVec4D? operator /(MVec4D? m, MVec4D? n) => Div(m, n); + public static MVec4D? operator /(double? c, MVec4D? n) => Div(c, n); + public static MVec4D? operator /(MVec4D? m, double? c) => Div(m, c); + + public static MVec4D? operator &(MVec4D? m, MVec4D? n) => Inner(m, n); + public static MVec4D? operator &(MVec4D? m, double? c) => Inner(m, c); + public static MVec4D? operator &(double? c, MVec4D? n) => Inner(c, n); + + public static MVec4D? operator ^(MVec4D? m, MVec4D? n) => Outer(m, n); + public static MVec4D? operator ^(MVec4D? m, double? c) => Outer(m, c); + public static MVec4D? operator ^(double? c, MVec4D? n) => Outer(c, n); + + + public static MVec4D operator -(MVec4D m) => m.Neg; + public static MVec4D? operator ~(MVec4D m) => m.Inv; + public static MVec4D operator !(MVec4D m) => m.Rev; + + public static implicit operator MVec4D(double v) => Vec0(v); + public static implicit operator MVec4D(Vector3d v) => Vec1(v.X, v.Y, v.Z, 0); + public static implicit operator MVec4D(Vector4d v) => Vec1(v.X, v.Y, v.Z, v.W); + #endregion #region Equality diff --git a/Tetrahedrons/TetrahedronWindow.cs b/Tetrahedrons/TetrahedronWindow.cs index d7061ed..07bec16 100644 --- a/Tetrahedrons/TetrahedronWindow.cs +++ b/Tetrahedrons/TetrahedronWindow.cs @@ -255,14 +255,14 @@ namespace Tetrahedrons if (Keyboard[Key.Left]) { - var rot = MVec3D.Complex(-e.Time, MVec3D.Unit12); + var rot = MVec3D.Rotor(-e.Time, MVec3D.Unit12); _b1 = rot * _b1 * ~rot; // todo make blade classes _b2 = rot * _b2 * ~rot; } if (Keyboard[Key.Right]) { - var rot = MVec3D.Complex(e.Time, MVec3D.Unit12); + var rot = MVec3D.Rotor(e.Time, MVec3D.Unit12); _b1 = rot * _b1 * ~rot; _b2 = rot * _b2 * ~rot; } From ac5cd8dd1dec9dd612e874134fd44446a939a491 Mon Sep 17 00:00:00 2001 From: allem Date: Tue, 3 Apr 2018 23:47:03 -0400 Subject: [PATCH 3/7] removed 3d code --- Tetrahedrons/MVec3d.cs | 195 ----------------- Tetrahedrons/Program.cs | 16 +- Tetrahedrons/TetrahedronWindow.cs | 353 ++++-------------------------- Tetrahedrons/Tetrahedrons.csproj | 1 - 4 files changed, 39 insertions(+), 526 deletions(-) delete mode 100644 Tetrahedrons/MVec3d.cs diff --git a/Tetrahedrons/MVec3d.cs b/Tetrahedrons/MVec3d.cs deleted file mode 100644 index e4c2cf9..0000000 --- a/Tetrahedrons/MVec3d.cs +++ /dev/null @@ -1,195 +0,0 @@ -using System; -using OpenTK; - -namespace Tetrahedrons -{ - public struct MVec3D - { - public static readonly MVec3D Zero = new MVec3D(0, 0, 0, 0, 0, 0, 0, 0); - public static readonly MVec3D One = new MVec3D(1, 1, 1, 1, 1, 1, 1, 1); - - public static readonly MVec3D Unit = new MVec3D(1, 0, 0, 0, 0, 0, 0, 0); - public static readonly MVec3D Unit1 = new MVec3D(0, 1, 0, 0, 0, 0, 0, 0); - public static readonly MVec3D Unit2 = new MVec3D(0, 0, 1, 0, 0, 0, 0, 0); - public static readonly MVec3D Unit3 = new MVec3D(0, 0, 0, 1, 0, 0, 0, 0); - public static readonly MVec3D Unit23 = new MVec3D(0, 0, 0, 0, 1, 0, 0, 0); - public static readonly MVec3D Unit31 = new MVec3D(0, 0, 0, 0, 0, 1, 0, 0); - public static readonly MVec3D Unit12 = new MVec3D(0, 0, 0, 0, 0, 0, 1, 0); - public static readonly MVec3D Unit123 = new MVec3D(0, 0, 0, 0, 0, 0, 0, 1); - - public double E; - public double E1; - public double E2; - public double E3; - public double E23; - public double E31; - public double E12; - public double E123; - - public MVec3D(double e, double e1, double e2, double e3, double e23, double e31, double e12, double e123) - { - E = e; - E1 = e1; - E2 = e2; - E3 = e3; - E23 = e23; - E31 = e31; - E12 = e12; - E123 = e123; - } - - public MVec3D(MVec3D s, MVec3D v, MVec3D b, MVec3D t) - : this(s.E, v.E1, v.E2, v.E3, b.E23, b.E31, b.E12, t.E123) - { - } - - public static MVec3D Scalar(double e) - { - return new MVec3D(e, 0, 0, 0, 0, 0, 0, 0); - } - - public static MVec3D Vector(double e1, double e2, double e3) - { - return new MVec3D(0, e1, e2, e3, 0, 0, 0, 0); - } - - public static MVec3D Vector(Vector3d v) - { - return new MVec3D(0, v.X, v.Y, v.Z, 0, 0, 0, 0); - } - - public static MVec3D Bivector(double e23, double e31, double e12) - { - return new MVec3D(0, 0, 0, 0, e23, e31, e12, 0); - } - - public static MVec3D Trivector(double e123) - { - return new MVec3D(0, 0, 0, 0, 0, 0, 0, e123); - } - - public static MVec3D Rotor(double angle, MVec3D plane) - { - return Math.Cos(angle) + Math.Sin(angle) * plane; - } - - public static MVec3D Add(MVec3D m, MVec3D n) - { - return new MVec3D(m.E + n.E, - m.E1 + n.E1, m.E2 + n.E2, m.E3 + n.E3, - m.E23 + n.E23, m.E31 + n.E31, m.E12 + n.E12, - m.E123 + n.E123); - } - - public static MVec3D Sub(MVec3D m, MVec3D n) - { - return new MVec3D(m.E - n.E, - m.E1 - n.E1, m.E2 - n.E2, m.E3 - n.E3, - m.E23 - n.E23, m.E31 - n.E31, m.E12 - n.E12, - m.E123 - n.E123); - } - - public static MVec3D Mult(MVec3D m, MVec3D n) - { - return new MVec3D( - m.E * n.E + m.E1 * n.E1 + m.E2 * n.E2 + m.E3 * n.E3 - m.E23 * n.E23 - m.E31 * n.E31 - m.E12 * n.E12 - m.E123 * n.E123, - m.E * n.E1 + m.E1 * n.E - m.E2 * n.E12 + m.E3 * n.E31 - m.E23 * n.E123 - m.E31 * n.E3 + m.E12 * n.E2 - m.E123 * n.E23, - m.E * n.E2 + m.E1 * n.E12 + m.E2 * n.E - m.E3 * n.E23 + m.E23 * n.E3 - m.E31 * n.E123 - m.E12 * n.E1 - m.E123 * n.E31, - m.E * n.E3 - m.E1 * n.E31 + m.E2 * n.E23 + m.E3 * n.E - m.E23 * n.E2 + m.E31 * n.E1 - m.E12 * n.E123 - m.E123 * n.E12, - m.E * n.E23 + m.E1 * n.E123 + m.E2 * n.E3 - m.E3 * n.E2 + m.E23 * n.E - m.E31 * n.E12 + m.E12 * n.E31 + m.E123 * n.E1, - m.E * n.E31 - m.E1 * n.E3 + m.E2 * n.E123 + m.E3 * n.E1 + m.E23 * n.E12 + m.E31 * n.E - m.E12 * n.E23 + m.E123 * n.E2, - m.E * n.E12 + m.E1 * n.E2 - m.E2 * n.E1 + m.E3 * n.E123 - m.E23 * n.E31 + m.E31 * n.E23 + m.E12 * n.E + m.E123 * n.E3, - m.E * n.E123 + m.E1 * n.E23 + m.E2 * n.E31 + m.E3 * n.E12 + m.E23 * n.E1 + m.E31 * n.E2 + m.E12 * n.E3 + m.E123 * n.E - ); - } - - public double Norm2 => E * E + E1 * E1 + E2 * E2 + E3 * E3 + E23 * E23 + E31 * E31 + E12 * E12 + E123 * E123; - public double Norm => Math.Sqrt(Norm2); - - public MVec3D Inv - { - get - { - var n2 = Norm2; - return new MVec3D(E / n2, -E1 / n2, -E2 / n2, -E3 / n2, -E23 / n2, -E31 / n2, -E12 / n2, E123 / n2); - } - } - - public MVec3D Grade(int k) - { - switch (k) - { - case 0: - return Scalar(E); - case 1: - return Vector(E1, E2, E3); - case 2: - return Bivector(E23, E31, E12); - case 3: - return Trivector(E123); - default: - return new MVec3D(); - } - } - - public static MVec3D Inner(MVec3D m, MVec3D n) - { - var m0 = m.Grade(0); - var m1 = m.Grade(1); - var m2 = m.Grade(2); - var m3 = m.Grade(3); - - var n0 = n.Grade(0); - var n1 = n.Grade(1); - var n2 = n.Grade(2); - var n3 = n.Grade(3); - - return new MVec3D( - m0 * n0 + m1 * n1 + m2 * n2 + m3 * n3, - m0 * n1 + m1 * n2 + m2 * n3, - m0 * n2 + m1 * n3, - m0 * n3 - ); - } - - public static MVec3D Outer(MVec3D m, MVec3D n) - { - var m0 = m.Grade(0); - var m1 = m.Grade(1); - var m2 = m.Grade(2); - var m3 = m.Grade(3); - - var n0 = n.Grade(0); - var n1 = n.Grade(1); - var n2 = n.Grade(2); - var n3 = n.Grade(3); - - return new MVec3D( - m0 * n0, - m1 * n0 + m0 * n1, - m2 * n0 + m1 * n1 + m0 * n2, - m3 * n0 + m2 * n1 + m1 * n2 + m0 * m3 - ); - } - - public static MVec3D operator ~(MVec3D m) => m.Inv; - public static MVec3D operator *(MVec3D m, MVec3D n) => Mult(m, n); - - public static MVec3D operator +(MVec3D m, MVec3D n) => Add(m, n); - public static MVec3D operator -(MVec3D m, MVec3D n) => Sub(m, n); - - public static MVec3D operator &(MVec3D m, MVec3D n) => Inner(m, n); - public static MVec3D operator ^(MVec3D m, MVec3D n) => Outer(m, n); - - public Vector3d VectorPart => new Vector3d(E1, E2, E3); - - public static implicit operator MVec3D(Vector3 v) => new MVec3D(0, v.X, v.Y, v.Z, 0, 0, 0, 0); - public static implicit operator MVec3D(Vector3d v) => new MVec3D(0, v.X, v.Y, v.Z, 0, 0, 0, 0); - public static implicit operator MVec3D(double s) => new MVec3D(s, 0, 0, 0, 0, 0, 0, 0); - - public override string ToString() - { - return $"({E:0.00} {E1:0.00}e1 {E2:0.00}e2 {E3:0.00}e3 {E23:0.00}e23 {E31:0.00}e31 {E12:0.00}e12 {E123:0.00}e123)"; - } - } -} \ No newline at end of file diff --git a/Tetrahedrons/Program.cs b/Tetrahedrons/Program.cs index 2015418..e2f16f6 100644 --- a/Tetrahedrons/Program.cs +++ b/Tetrahedrons/Program.cs @@ -5,24 +5,10 @@ namespace Tetrahedrons { internal class Program : GameWindow { - public static void NotMain(string[] args) + public static void Main(string[] args) { using (var p = new TetrahedronWindow()) p.Run(); } - - public static void Main(string[] args) - { - var m = new MVec4D(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0); - var i = m.Inv ?? MVec4D.Zero; - - Console.Out.WriteLine("m = {0}", m); - Console.Out.WriteLine("m.Rev = {0}", m.Rev); - Console.Out.WriteLine("MVec4D.Mul(m, m.Rev) = {0}", MVec4D.Mul(m, m.Rev)); - Console.Out.WriteLine("m.Inv = {0}", i); - - Console.Out.WriteLine("MVec4D.Mul(m, m.Inv) = {0}", MVec4D.Mul(m, i)); - Console.Out.WriteLine("MVec4D.Mul(m.Inv, m) = {0}", MVec4D.Mul(i, m)); - } } } \ No newline at end of file diff --git a/Tetrahedrons/TetrahedronWindow.cs b/Tetrahedrons/TetrahedronWindow.cs index 07bec16..0e26716 100644 --- a/Tetrahedrons/TetrahedronWindow.cs +++ b/Tetrahedrons/TetrahedronWindow.cs @@ -10,332 +10,55 @@ using OpenTK.Input; namespace Tetrahedrons { - public struct Polygon - { - public Vector3d[] Points; - public Color Color; + public class TetrahedronWindow : GameWindow + { + private Matrix4 _proj3d; + private Matrix4 _view; - public Polygon(Vector3d[] points) : this(points, Color.DodgerBlue) - { - } + Vector3d a = Vector3d.UnitX; + Vector3d b = Vector3d.UnitY; - public Polygon(Vector3d[] points, Color color) - { - Points = points; - Color = color; - } + protected override void OnLoad(EventArgs e) + { + base.OnLoad(e); - public void Draw() - { - GL.Begin(PrimitiveType.TriangleStrip); - GL.Color3(Color); - foreach (var v in Points) - GL.Vertex3(v); - GL.End(); - } + X = (DisplayDevice.Default.Width - Width) / 2; + Y = (DisplayDevice.Default.Height - Height) / 2; - public static Polygon operator +(Polygon p, MVec3D m) - { - var pts = new Vector3d[p.Points.Length]; - for (var i = 0; i < p.Points.Length; i++) - pts[i] = (p.Points[i] + m).VectorPart; - return new Polygon(pts, p.Color); - } + _view = Matrix4.LookAt(Vector3.Zero, -new Vector3(1, .6f, 1f), Vector3.UnitZ); + } - public static Polygon operator -(Polygon p, MVec3D m) - { - var pts = new Vector3d[p.Points.Length]; - for (var i = 0; i < p.Points.Length; i++) - pts[i] = (p.Points[i] - m).VectorPart; - return new Polygon(pts, p.Color); - } + protected override void OnRenderFrame(FrameEventArgs e) + { + base.OnRenderFrame(e); - public static Polygon operator *(Polygon p, MVec3D m) - { - var pts = new Vector3d[p.Points.Length]; - for (var i = 0; i < p.Points.Length; i++) - pts[i] = (p.Points[i] * m).VectorPart; - return new Polygon(pts, p.Color); - } + GL.Viewport(ClientRectangle); - public static Polygon operator ^(Polygon p, MVec3D m) - { - var pts = new Vector3d[p.Points.Length]; - for (var i = 0; i < p.Points.Length; i++) - pts[i] = (p.Points[i] ^ m).VectorPart; - return new Polygon(pts, p.Color); - } + GL.Clear(ClearBufferMask.ColorBufferBit | ClearBufferMask.DepthBufferBit); + GL.PointSize(10f); - public static Polygon operator &(Polygon p, MVec3D m) - { - var pts = new Vector3d[p.Points.Length]; - for (var i = 0; i < p.Points.Length; i++) - pts[i] = (p.Points[i] & m).VectorPart; - return new Polygon(pts, p.Color); - } - } + GL.PushMatrix(); - public class Tetrahedron - { - public readonly Vector3d[] Points; - public Color Color; - public Color WireColor; + GL.MatrixMode(MatrixMode.Projection); + GL.LoadMatrix(ref _proj3d); + GL.MatrixMode(MatrixMode.Modelview); + GL.LoadMatrix(ref _view); - public Tetrahedron(Vector3d p1, Vector3d p2, Vector3d p3, Vector3d p4) : this(p1, p2, p3, p4, Color.GhostWhite) - { - } + GL.Begin(PrimitiveType.TriangleStrip); + GL.Vertex3(-a - b); + GL.Vertex3(-a + b); + GL.Vertex3(a - b); + GL.Vertex3(a + b); + GL.End(); - public Tetrahedron(Vector3d p1, Vector3d p2, Vector3d p3, Vector3d p4, Color color) : this(p1, p2, p3, p4, color, Color.Maroon) - { - } + SwapBuffers(); + } - public Tetrahedron(Vector3d p1, Vector3d p2, Vector3d p3, Vector3d p4, Color color, Color wireColor) - { - Color = color; - WireColor = wireColor; - Points = new[] {p1, p2, p3, p4}; - } + protected override void OnUpdateFrame(FrameEventArgs e) + { + base.OnUpdateFrame(e); - public void Draw(bool wire = true) - { - if (wire) - { - GL.Begin(PrimitiveType.Lines); - GL.Color3(WireColor); - - for (var i = 0; i < 4; i++) - for (var j = i + 1; j < 4; j++) - { - GL.Vertex3(Points[i]); - GL.Vertex3(Points[j]); - } - - GL.End(); - } - else - { - GL.Begin(PrimitiveType.Triangles); - GL.Color3(Color); - - for (var i = 0; i < 4; i++) - for (var j = i + 1; j < 4; j++) - for (var k = j + 1; k < 4; k++) - { - GL.Vertex3(Points[i]); - GL.Vertex3(Points[j]); - GL.Vertex3(Points[k]); - } - - GL.End(); - } - } - - public Polygon Intersection(MVec3D blade, Vector3d pivot) - { - var pts = new Vector3d[4]; - for (var i = 0; i < 4; i++) - { - pts[i] = Points[i] - pivot; - } - - var a = blade * MVec3D.Unit123; - var sides = new double[4]; - for (var i = 0; i < 4; i++) - { - sides[i] = (a & ((pts[i] ^ blade) * ~blade)).E; // normal dot rejection - } - - var vecs = new List(4); - - for (var i = 0; i < 4; i++) - { - for (var j = i + 1; j < 4; j++) - { - if (sides[i] * sides[j] >= 0) continue; // both pts on same side - - var t = ((blade ^ pts[i]) * ~(blade ^ (pts[j] - pts[i]))).E; - var point = pts[i] + t * (pts[j] - pts[i]); - vecs.Add(point + pivot); - } - } - - return new Polygon(vecs.ToArray()); - } - } - - public class TetrahedronWindow : GameWindow - { - private Matrix4 _proj3d; - private Matrix4 _proj2d; - private Matrix4 _view; - - private double _t; - private bool _pause; - private bool _div; - private Tetrahedron _tetra; - private Polygon _polyg; - - private MVec3D _b1 = MVec3D.Unit1; - private MVec3D _b2 = MVec3D.Unit3; - private MVec3D _pivot = MVec3D.Zero; - - private Polygon _square = new Polygon(new[] {new Vector3d(-1, -1, 0), new Vector3d(-1, 1, 0), new Vector3d(1, -1, 0), new Vector3d(1, 1, 0),}, Color.Red); - - protected override void OnLoad(EventArgs e) - { - base.OnLoad(e); - - X = (DisplayDevice.Default.Width - Width) / 2; - Y = (DisplayDevice.Default.Height - Height) / 2; - - _tetra = new Tetrahedron(new Vector3d(1, 1, 1), new Vector3d(-1, -1, 1), new Vector3d(-1, 1, -1), new Vector3d(1, -1, -1)); - _polyg = new Polygon(new[] {new Vector3d(-1, -1, 0), new Vector3d(-1, 1, 0), new Vector3d(1, -1, 0), new Vector3d(1, 1, 0)}); - - _view = Matrix4.LookAt(Vector3.Zero, -new Vector3(1, .6f, 1f), Vector3.UnitZ); - } - - protected override void OnRenderFrame(FrameEventArgs e) - { - base.OnRenderFrame(e); - - GL.Viewport(ClientRectangle); - - GL.Clear(ClearBufferMask.ColorBufferBit | ClearBufferMask.DepthBufferBit); - GL.PointSize(10f); - - GL.PushMatrix(); - - if (_div) - { - GL.MatrixMode(MatrixMode.Projection); - GL.LoadMatrix(ref _proj2d); - - _polyg.Draw(); - } - else - { - GL.MatrixMode(MatrixMode.Projection); - GL.LoadMatrix(ref _proj3d); - GL.MatrixMode(MatrixMode.Modelview); - GL.LoadMatrix(ref _view); - - - GL.Enable(EnableCap.DepthTest); - _tetra.Draw(wire: false); - GL.Disable(EnableCap.DepthTest); - _polyg.Draw(); - - GL.Begin(PrimitiveType.Points); - GL.Color3(Color.DarkBlue); - GL.Vertex3(_pivot.VectorPart); - GL.End(); - - GL.Begin(PrimitiveType.Lines); - GL.Vertex3(_pivot.VectorPart); - GL.Vertex3((_pivot + (_b2 ^ _b1) * MVec3D.Unit123 * .5).VectorPart); - GL.End(); - - GL.Enable(EnableCap.DepthTest); - _tetra.Draw(wire: true); - } - - GL.PopMatrix(); - - SwapBuffers(); - } - - protected override void OnUpdateFrame(FrameEventArgs e) - { - base.OnUpdateFrame(e); - - _proj3d = Matrix4.CreateOrthographic(6, 6f * Height / Width, -2, 2); - _proj2d = Matrix4.CreateOrthographic(4, 4f * Height / Width, -1, 1); - - _t += e.Time; - - if (Keyboard[Key.Left]) - { - var rot = MVec3D.Rotor(-e.Time, MVec3D.Unit12); - _b1 = rot * _b1 * ~rot; // todo make blade classes - _b2 = rot * _b2 * ~rot; - } - - if (Keyboard[Key.Right]) - { - var rot = MVec3D.Rotor(e.Time, MVec3D.Unit12); - _b1 = rot * _b1 * ~rot; - _b2 = rot * _b2 * ~rot; - } - - if (Keyboard[Key.Up]) - { - _pivot -= (_b1 ^ _b2) * MVec3D.Unit123 * e.Time; - } - - if (Keyboard[Key.Down]) - { - _pivot += (_b1 ^ _b2) * MVec3D.Unit123 * e.Time; - } - - if (_pause) - return; - -// var b1 = MVec3d.Unit2; -// var b2 = MVec3d.Unit3; -// var pivot = new Vector3d(Math.Cos(_t), 0, 0); - -// var b1 = MVec3d.Vector(0, Math.Cos(_t * .2), Math.Sin(_t * .5)); -// var b2 = MVec3d.Vector(Math.Cos(_t * .6), 0, Math.Sin(_t * .4)); -// var pivot = Vector3d.Zero; - -// var b1 = MVec3d.Vector(Math.Cos(_t/2), Math.Sin(_t/2), 0); -// var b2 = MVec3d.Unit3; -// var pivot = Vector3d.Zero; - -// var blade = MVec3d.Outer(b1, b2); - -// var rot = MVec3d.Complex(_t / 2, MVec3d.Unit12); -// var b1 = rot * MVec3d.Unit2 * ~rot; -// var b2 = rot * MVec3d.Unit3 * ~rot; -// var pln = MVec3d.Outer(b1, b2); -// var blade = pln; -// var pivot = Vector3d.Zero; - - _polyg = _tetra.Intersection(_b1 ^ _b2, _pivot.VectorPart); - - if (_div) - { - var px = ((_pivot & _b1) * ~_b1) & _b1; - var py = ((_pivot ^ _b1) * ~_b1) & ((_b2 ^ _b1) * ~_b1); - - for (var i = 0; i < _polyg.Points.Length; i++) - { - MVec3D p = _polyg.Points[i]; - var x = ((p & _b1) * ~_b1) & _b1; - var y = ((p ^ _b1) * ~_b1) & ((_b2 ^ _b1) * ~_b1); - _polyg.Points[i] = new Vector3d(x.E - px.E, y.E - py.E, 0); - } - } - } - - protected override void OnKeyDown(KeyboardKeyEventArgs e) - { - base.OnKeyDown(e); - - switch (e.Key) - { - case Key.Space: - _pause = !_pause; - break; - case Key.A: - _div = !_div; - break; - case Key.R: - _b1 = MVec3D.Unit1; - _b2 = MVec3D.Unit3; - _pivot = MVec3D.Zero; - break; - } - } - } + _proj3d = Matrix4.CreateOrthographic(6, 6f * Height / Width, -2, 2); + } + } } \ No newline at end of file diff --git a/Tetrahedrons/Tetrahedrons.csproj b/Tetrahedrons/Tetrahedrons.csproj index 70b1a79..ad503ef 100644 --- a/Tetrahedrons/Tetrahedrons.csproj +++ b/Tetrahedrons/Tetrahedrons.csproj @@ -42,7 +42,6 @@ - From 632172e8769b4ffba9b9c14547e62b68966ccccd Mon Sep 17 00:00:00 2001 From: allem Date: Wed, 4 Apr 2018 02:04:11 -0400 Subject: [PATCH 4/7] simple projections through 4d --- Tetrahedrons/MVec4D.cs | 3 ++ Tetrahedrons/TetrahedronWindow.cs | 56 +++++++++++++++++++++++++++---- 2 files changed, 52 insertions(+), 7 deletions(-) diff --git a/Tetrahedrons/MVec4D.cs b/Tetrahedrons/MVec4D.cs index ddc59be..eb9bdd5 100644 --- a/Tetrahedrons/MVec4D.cs +++ b/Tetrahedrons/MVec4D.cs @@ -61,6 +61,9 @@ namespace Tetrahedrons } } + public Vector4d V4 => new Vector4d(X, Y, Z, W); + public Vector3d V3 => new Vector3d(X, Y, Z); + #region Swizzle public double Yx diff --git a/Tetrahedrons/TetrahedronWindow.cs b/Tetrahedrons/TetrahedronWindow.cs index 0e26716..b95f896 100644 --- a/Tetrahedrons/TetrahedronWindow.cs +++ b/Tetrahedrons/TetrahedronWindow.cs @@ -1,6 +1,7 @@ using System; using System.Collections.Generic; using System.Drawing; +using System.Drawing.Printing; using System.Net.Mail; using System.Security; using OpenTK; @@ -15,8 +16,10 @@ namespace Tetrahedrons private Matrix4 _proj3d; private Matrix4 _view; - Vector3d a = Vector3d.UnitX; - Vector3d b = Vector3d.UnitY; + MVec4D a = MVec4D.UnitX; + MVec4D b = MVec4D.UnitY; + MVec4D v = MVec4D.Vec1(1, 1, 0, -1); + MVec4D t = MVec4D.Vec1(0, 0, 0, 1); protected override void OnLoad(EventArgs e) { @@ -25,7 +28,8 @@ namespace Tetrahedrons X = (DisplayDevice.Default.Width - Width) / 2; Y = (DisplayDevice.Default.Height - Height) / 2; - _view = Matrix4.LookAt(Vector3.Zero, -new Vector3(1, .6f, 1f), Vector3.UnitZ); + _view = Matrix4.LookAt(Vector3.Zero, -new Vector3(1, 1, 1), Vector3.UnitZ); + } protected override void OnRenderFrame(FrameEventArgs e) @@ -35,7 +39,9 @@ namespace Tetrahedrons GL.Viewport(ClientRectangle); GL.Clear(ClearBufferMask.ColorBufferBit | ClearBufferMask.DepthBufferBit); + GL.Enable(EnableCap.DepthTest); GL.PointSize(10f); + GL.LineWidth(3f); GL.PushMatrix(); @@ -44,11 +50,41 @@ namespace Tetrahedrons GL.MatrixMode(MatrixMode.Modelview); GL.LoadMatrix(ref _view); + GL.Disable(EnableCap.DepthTest); GL.Begin(PrimitiveType.TriangleStrip); - GL.Vertex3(-a - b); - GL.Vertex3(-a + b); - GL.Vertex3(a - b); - GL.Vertex3(a + b); + GL.Color3(1d, 1, 1); + GL.Vertex3((-a - b).V3); + GL.Vertex3((-a + b).V3); + GL.Vertex3((a - b).V3); + GL.Vertex3((a + b).V3); + GL.End(); + GL.Enable(EnableCap.DepthTest); + + GL.Begin(PrimitiveType.Lines); + GL.Color3(1d,0,0); GL.Vertex3(0,0,0); GL.Vertex3(1,0,0); + GL.Color3(0,1d,0); GL.Vertex3(0,0,0); GL.Vertex3(0,1,0); + GL.Color3(0,0,1d); GL.Vertex3(0,0,0); GL.Vertex3(0,0,1); + GL.End(); + + var B = MVec4D.UnitXyz; + var bt = B ^ t; + var bv = B ^ v; + var alph = bt / bv ?? MVec4D.Zero; + var dt = alph * v; + var p = t - dt; + + Console.Out.WriteLine("p = {0}", p); + + GL.Begin(PrimitiveType.Lines); + GL.Color3(.5 + t.W/4, 0, .5 - t.W/4); + GL.Vertex3(t.V3); + GL.Color3(.5 + (t + v).W/4, 0, .5 - (t + v).W/4); + GL.Vertex3((t + v).V3); + GL.End(); + + GL.Begin(PrimitiveType.Points); + GL.Color3(0d, 0, 1); + GL.Vertex3(p.V3); GL.End(); SwapBuffers(); @@ -59,6 +95,12 @@ namespace Tetrahedrons base.OnUpdateFrame(e); _proj3d = Matrix4.CreateOrthographic(6, 6f * Height / Width, -2, 2); + + var rv = MVec4D.Rotor(e.Time, MVec4D.UnitXy); + v = rv * v * !rv; + + var rt = MVec4D.Rotor(e.Time / 4, MVec4D.UnitZw); + t = rt * t * !rt; } } } \ No newline at end of file From 06096bf1aadafccd1ed3d3a713df79696bb2d1ba Mon Sep 17 00:00:00 2001 From: allem Date: Wed, 4 Apr 2018 14:06:57 -0400 Subject: [PATCH 5/7] buggy 4-simplex intersections --- Tetrahedrons/MVec4D.cs | 15 +---- Tetrahedrons/Simplex.cs | 90 ++++++++++++++++++++++++++++ Tetrahedrons/TetrahedronWindow.cs | 98 +++++++++++++++++++------------ Tetrahedrons/Tetrahedrons.csproj | 2 + Tetrahedrons/Util.cs | 25 ++++++++ 5 files changed, 181 insertions(+), 49 deletions(-) create mode 100644 Tetrahedrons/Simplex.cs create mode 100644 Tetrahedrons/Util.cs diff --git a/Tetrahedrons/MVec4D.cs b/Tetrahedrons/MVec4D.cs index eb9bdd5..bfb056b 100644 --- a/Tetrahedrons/MVec4D.cs +++ b/Tetrahedrons/MVec4D.cs @@ -2,6 +2,7 @@ using System.Collections.Generic; using System.Data.SqlClient; using System.Dynamic; +using System.Runtime.CompilerServices; using System.Runtime.Remoting.Messaging; using System.Text; using OpenTK; @@ -42,6 +43,8 @@ namespace Tetrahedrons public double Norm => Math.Sqrt(Norm2); + public MVec4D Normalized => this * (1 / Norm); + public MVec4D Grade(int k) { switch (k) @@ -929,12 +932,8 @@ namespace Tetrahedrons ); } - public static MVec4D? Div(MVec4D m, MVec4D n) => Mul(m, n.Inv); - public static MVec4D Div(MVec4D m, double c) => Mul(m, 1 / c); - public static MVec4D? Div(double c, MVec4D n) => Mul(c, n.Inv); - public static MVec4D Inner(MVec4D m, MVec4D n) { var m0 = m.Grade(0); @@ -1030,10 +1029,6 @@ namespace Tetrahedrons public static MVec4D operator *(double c, MVec4D n) => Mul(c, n); public static MVec4D operator *(MVec4D m, double c) => Mul(m, c); - public static MVec4D? operator /(MVec4D m, MVec4D n) => Div(m, n); - public static MVec4D? operator /(double c, MVec4D n) => Div(c, n); - public static MVec4D operator /(MVec4D m, double c) => Div(m, c); - public static MVec4D operator &(MVec4D m, MVec4D n) => Inner(m, n); public static MVec4D operator &(MVec4D m, double c) => Inner(m, c); public static MVec4D operator &(double c, MVec4D n) => Inner(c, n); @@ -1054,10 +1049,6 @@ namespace Tetrahedrons public static MVec4D? operator *(double? c, MVec4D? n) => Mul(c, n); public static MVec4D? operator *(MVec4D? m, double? c) => Mul(m, c); - public static MVec4D? operator /(MVec4D? m, MVec4D? n) => Div(m, n); - public static MVec4D? operator /(double? c, MVec4D? n) => Div(c, n); - public static MVec4D? operator /(MVec4D? m, double? c) => Div(m, c); - public static MVec4D? operator &(MVec4D? m, MVec4D? n) => Inner(m, n); public static MVec4D? operator &(MVec4D? m, double? c) => Inner(m, c); public static MVec4D? operator &(double? c, MVec4D? n) => Inner(c, n); diff --git a/Tetrahedrons/Simplex.cs b/Tetrahedrons/Simplex.cs new file mode 100644 index 0000000..1fc6e42 --- /dev/null +++ b/Tetrahedrons/Simplex.cs @@ -0,0 +1,90 @@ +using System.Collections.Generic; +using System.Linq; + +namespace Tetrahedrons +{ + public class Simplex + { + public MVec4D[] Verts; + + public Simplex(IEnumerable verts) + { + Verts = verts.ToArray(); + } + + public Simplex(params MVec4D[] verts) + { + Verts = verts; + } + + public bool[] Sides(MVec4D blade, MVec4D pivot) + { + var sides = new bool[Verts.Length]; + + var perp = blade.Dual; + + for (var i = 0; i < Verts.Length; i++) + { + var m = Verts[i]; + var p = m - pivot; + var rej = (p ^ blade) * ~blade; + if (!rej.HasValue) continue; + var dot = rej & perp; + + sides[i] = dot.Value.S > 0; + } + + return sides; + } + + public Simplex Intersect(MVec4D blade, MVec4D pivot) + { + var sides = Sides(blade, pivot); + var verts = new List(); + + for (var i = 0; i < Verts.Length; i++) + for (var j = i + 1; j < Verts.Length; j++) + { + if (sides[i] == sides[j]) continue; + var m = Verts[i] - pivot; + var n = Verts[j] - pivot; + + var v = m - n; + var alph = (blade ^ m) * ~(blade ^ v); + if (!alph.HasValue) continue; + var p = m - v * alph.Value; + verts.Add(p); + } + + return new Simplex(verts); + } + + public int[] Faces() + { + switch (Verts.Length) + { + case 3: + return new[] {0, 1, 2}; + case 4: + return new[] + { + 0, 1, 2, + 0, 1, 3, + 0, 2, 3, + 1, 2, 3 + }; + case 6: + return new[] + { + 0, 2, 4 + // todo: vertices of a face must be on lines which share a vertex in the simplex + // that's the pattern, and why it "just works" for the 3-simplex. + // think of it as truncating the embedded simplex + }; + + default: + return new int[0]; + } + } + } +} \ No newline at end of file diff --git a/Tetrahedrons/TetrahedronWindow.cs b/Tetrahedrons/TetrahedronWindow.cs index b95f896..124e9da 100644 --- a/Tetrahedrons/TetrahedronWindow.cs +++ b/Tetrahedrons/TetrahedronWindow.cs @@ -1,13 +1,17 @@ using System; -using System.Collections.Generic; +using System.Diagnostics.PerformanceData; using System.Drawing; using System.Drawing.Printing; +using System.Net; using System.Net.Mail; +using System.Runtime.InteropServices.ComTypes; using System.Security; +using System.Xml.XPath; using OpenTK; using OpenTK.Audio.OpenAL; using OpenTK.Graphics.OpenGL; using OpenTK.Input; +using OpenTK.Platform.Windows; namespace Tetrahedrons { @@ -16,10 +20,13 @@ namespace Tetrahedrons private Matrix4 _proj3d; private Matrix4 _view; - MVec4D a = MVec4D.UnitX; - MVec4D b = MVec4D.UnitY; - MVec4D v = MVec4D.Vec1(1, 1, 0, -1); - MVec4D t = MVec4D.Vec1(0, 0, 0, 1); + private MVec4D _a = MVec4D.UnitX; + private MVec4D _b = MVec4D.UnitY; + + private Simplex _pent; + private Simplex _intr; + + private double _t; protected override void OnLoad(EventArgs e) { @@ -28,8 +35,14 @@ namespace Tetrahedrons X = (DisplayDevice.Default.Width - Width) / 2; Y = (DisplayDevice.Default.Height - Height) / 2; - _view = Matrix4.LookAt(Vector3.Zero, -new Vector3(1, 1, 1), Vector3.UnitZ); + _view = Matrix4.LookAt(Vector3.Zero, -new Vector3(.86f, .5f, 1), Vector3.UnitZ); + _pent = new Simplex( // not-quite-regular pentatope + new Vector4d(1, -1, -1, -1), + new Vector4d(-1, 1, -1, -1), + new Vector4d(-1, -1, 1, -1), + new Vector4d(-1, -1, -1, 1), + new Vector4d(1, 1, 1, 1)); } protected override void OnRenderFrame(FrameEventArgs e) @@ -39,7 +52,6 @@ namespace Tetrahedrons GL.Viewport(ClientRectangle); GL.Clear(ClearBufferMask.ColorBufferBit | ClearBufferMask.DepthBufferBit); - GL.Enable(EnableCap.DepthTest); GL.PointSize(10f); GL.LineWidth(3f); @@ -50,41 +62,42 @@ namespace Tetrahedrons GL.MatrixMode(MatrixMode.Modelview); GL.LoadMatrix(ref _view); - GL.Disable(EnableCap.DepthTest); GL.Begin(PrimitiveType.TriangleStrip); GL.Color3(1d, 1, 1); - GL.Vertex3((-a - b).V3); - GL.Vertex3((-a + b).V3); - GL.Vertex3((a - b).V3); - GL.Vertex3((a + b).V3); + Util.Vertex3(-_a - _b); + Util.Vertex3(-_a + _b); + Util.Vertex3(_a - _b); + Util.Vertex3(_a + _b); GL.End(); + GL.Enable(EnableCap.DepthTest); - - GL.Begin(PrimitiveType.Lines); - GL.Color3(1d,0,0); GL.Vertex3(0,0,0); GL.Vertex3(1,0,0); - GL.Color3(0,1d,0); GL.Vertex3(0,0,0); GL.Vertex3(0,1,0); - GL.Color3(0,0,1d); GL.Vertex3(0,0,0); GL.Vertex3(0,0,1); + GL.Begin(PrimitiveType.Triangles); + foreach (var f in _intr.Faces()) + { + Util.Color3(_intr.Verts[f]); + Util.Vertex3(_intr.Verts[f]); + } + GL.End(); + GL.Disable(EnableCap.DepthTest); + + GL.Begin(PrimitiveType.Points); + foreach (var m in _pent.Verts) + Util.Vertex4(m); GL.End(); - var B = MVec4D.UnitXyz; - var bt = B ^ t; - var bv = B ^ v; - var alph = bt / bv ?? MVec4D.Zero; - var dt = alph * v; - var p = t - dt; - - Console.Out.WriteLine("p = {0}", p); - GL.Begin(PrimitiveType.Lines); - GL.Color3(.5 + t.W/4, 0, .5 - t.W/4); - GL.Vertex3(t.V3); - GL.Color3(.5 + (t + v).W/4, 0, .5 - (t + v).W/4); - GL.Vertex3((t + v).V3); + for (var i = 0; i < _pent.Verts.Length; i++) + for (var j = i + 1; j < _pent.Verts.Length; j++) + { + Util.Vertex4(_pent.Verts[i]); + Util.Vertex4(_pent.Verts[j]); + } + GL.End(); GL.Begin(PrimitiveType.Points); - GL.Color3(0d, 0, 1); - GL.Vertex3(p.V3); + foreach (var m in _intr.Verts) + Util.Vertex4(m); GL.End(); SwapBuffers(); @@ -94,13 +107,24 @@ namespace Tetrahedrons { base.OnUpdateFrame(e); - _proj3d = Matrix4.CreateOrthographic(6, 6f * Height / Width, -2, 2); + _proj3d = Matrix4.CreateOrthographic(6, 6f * Height / Width, -4, 4); - var rv = MVec4D.Rotor(e.Time, MVec4D.UnitXy); - v = rv * v * !rv; + _t += e.Time; - var rt = MVec4D.Rotor(e.Time / 4, MVec4D.UnitZw); - t = rt * t * !rt; + var pln = MVec4D.UnitXy + MVec4D.UnitZw; + + var r = MVec4D.Rotor(e.Time / 10, pln.Normalized); + for (var i = 0; i < _pent.Verts.Length; i++) + { + var m = _pent.Verts[i]; + m = r * m * !r; + _pent.Verts[i] = m; + } + + var blade = MVec4D.UnitXyz; +// var pivot = .5 * MVec4D.UnitW; + var pivot = 1.1 * Math.Sin(_t) * MVec4D.UnitW; + _intr = _pent.Intersect(blade, pivot); } } } \ No newline at end of file diff --git a/Tetrahedrons/Tetrahedrons.csproj b/Tetrahedrons/Tetrahedrons.csproj index ad503ef..546af1a 100644 --- a/Tetrahedrons/Tetrahedrons.csproj +++ b/Tetrahedrons/Tetrahedrons.csproj @@ -45,7 +45,9 @@ + + diff --git a/Tetrahedrons/Util.cs b/Tetrahedrons/Util.cs new file mode 100644 index 0000000..6f3aa54 --- /dev/null +++ b/Tetrahedrons/Util.cs @@ -0,0 +1,25 @@ +using OpenTK.Graphics.OpenGL; + +namespace Tetrahedrons +{ + public static class Util + { + public static void Color3(MVec4D m) + { + m *= .5; + GL.Color3(.5 + m.X, .5 + m.Y, .5 * m.Z); + } + + public static void Vertex4(MVec4D m) + { + var i = m.W / 2; + GL.Color3(.5 + i, 0, .5 - i); + GL.Vertex3(m.X, m.Y, m.Z); + } + + public static void Vertex3(MVec4D m) + { + GL.Vertex3(m.X, m.Y, m.Z); + } + } +} \ No newline at end of file From cb90af4fb71e4ff95c959cd4492b3a7608d587e0 Mon Sep 17 00:00:00 2001 From: allem Date: Wed, 4 Apr 2018 15:24:16 -0400 Subject: [PATCH 6/7] implemented 4-simplex intersections --- Tetrahedrons/MVec4D.cs | 10 ++- Tetrahedrons/Simplex.cs | 136 +++++++++++++++++++++--------- Tetrahedrons/TetrahedronWindow.cs | 60 ++++++------- 3 files changed, 133 insertions(+), 73 deletions(-) diff --git a/Tetrahedrons/MVec4D.cs b/Tetrahedrons/MVec4D.cs index bfb056b..9a0d8aa 100644 --- a/Tetrahedrons/MVec4D.cs +++ b/Tetrahedrons/MVec4D.cs @@ -988,7 +988,8 @@ namespace Tetrahedrons public static MVec4D Outer(double c, MVec4D n) => Mul(c, n); - + public static MVec4D Transform(MVec4D m, MVec4D t) => Mul(Mul(t, m), t.Rev); + public static MVec4D? Add(MVec4D? m, MVec4D? n) => m.HasValue && n.HasValue ? (MVec4D?) Add(m.Value, n.Value) : null; public static MVec4D? Add(MVec4D? m, double? n) => m.HasValue && n.HasValue ? (MVec4D?) Add(m.Value, n.Value) : null; public static MVec4D? Add(double? m, MVec4D? n) => n.HasValue && m.HasValue ? (MVec4D?) Add(m.Value, n.Value) : null; @@ -1013,6 +1014,8 @@ namespace Tetrahedrons public static MVec4D? Outer(MVec4D? m, double? n) => m.HasValue && n.HasValue ? (MVec4D?) Outer(m.Value, n.Value) : null; public static MVec4D? Outer(double? m, MVec4D? n) => n.HasValue && m.HasValue ? (MVec4D?) Outer(m.Value, n.Value) : null; + public static MVec4D? Transform(MVec4D? m, MVec4D? t) => Mul(Mul(t, m), t?.Rev); + #endregion #region Operators @@ -1037,6 +1040,9 @@ namespace Tetrahedrons public static MVec4D operator ^(MVec4D m, double c) => Outer(m, c); public static MVec4D operator ^(double c, MVec4D n) => Outer(c, n); + public static MVec4D operator |(MVec4D m, MVec4D t) => Transform(m, t); + + public static MVec4D? operator +(MVec4D? m, MVec4D? n) => Add(m, n); public static MVec4D? operator +(double? c, MVec4D? n) => Add(c, n); public static MVec4D? operator +(MVec4D? m, double? c) => Add(m, c); @@ -1057,7 +1063,9 @@ namespace Tetrahedrons public static MVec4D? operator ^(MVec4D? m, double? c) => Outer(m, c); public static MVec4D? operator ^(double? c, MVec4D? n) => Outer(c, n); + public static MVec4D? operator |(MVec4D? m, MVec4D? t) => Transform(m, t); + public static MVec4D operator -(MVec4D m) => m.Neg; public static MVec4D? operator ~(MVec4D m) => m.Inv; public static MVec4D operator !(MVec4D m) => m.Rev; diff --git a/Tetrahedrons/Simplex.cs b/Tetrahedrons/Simplex.cs index 1fc6e42..13d9cc1 100644 --- a/Tetrahedrons/Simplex.cs +++ b/Tetrahedrons/Simplex.cs @@ -1,5 +1,7 @@ -using System.Collections.Generic; +using System; +using System.Collections.Generic; using System.Linq; +using System.Net.Mail; namespace Tetrahedrons { @@ -7,14 +9,42 @@ namespace Tetrahedrons { public MVec4D[] Verts; + public int[] Faces; + public int[] Edjes; + public Simplex(IEnumerable verts) + : this(verts.ToArray()) { - Verts = verts.ToArray(); } public Simplex(params MVec4D[] verts) { - Verts = verts; + Verts = verts.ToArray(); + + var e = new List(); + for (var i = 0; i < Verts.Length; i++) + for (var j = i + 1; j < Verts.Length; j++) + { + e.Add(i); + e.Add(j); + } + + Edjes = e.ToArray(); + + Faces = new int[0]; + } + + public Simplex(IEnumerable verts, IEnumerable faces) + : this(verts) + { + Faces = faces.ToArray(); + } + + public Simplex(IEnumerable verts, IEnumerable faces, IEnumerable edges) + { + Verts = verts.ToArray(); + Faces = faces.ToArray(); + Edjes = edges.ToArray(); } public bool[] Sides(MVec4D blade, MVec4D pivot) @@ -41,50 +71,78 @@ namespace Tetrahedrons { var sides = Sides(blade, pivot); var verts = new List(); + var edges = new List(); + + var face_verts = new List[Verts.Length]; + for (var i = 0; i < face_verts.Length; i++) + face_verts[i] = new List(); - for (var i = 0; i < Verts.Length; i++) - for (var j = i + 1; j < Verts.Length; j++) { - if (sides[i] == sides[j]) continue; - var m = Verts[i] - pivot; - var n = Verts[j] - pivot; + var k = 0; + for (var i = 0; i < Verts.Length; i++) + for (var j = i + 1; j < Verts.Length; j++) + { + if (sides[i] == sides[j]) continue; + var m = Verts[i] - pivot; + var n = Verts[j] - pivot; - var v = m - n; - var alph = (blade ^ m) * ~(blade ^ v); - if (!alph.HasValue) continue; - var p = m - v * alph.Value; - verts.Add(p); + var v = m - n; + var alph = (blade ^ m) * ~(blade ^ v); + if (!alph.HasValue) continue; + var p = m - v * alph.Value; + + verts.Add(p); + face_verts[i].Add(k); + face_verts[j].Add(k); + + k++; + } } - return new Simplex(verts); - } - - public int[] Faces() - { - switch (Verts.Length) + var faces = new List(); + foreach (var vlst in face_verts) { - case 3: - return new[] {0, 1, 2}; - case 4: - return new[] - { - 0, 1, 2, - 0, 1, 3, - 0, 2, 3, - 1, 2, 3 - }; - case 6: - return new[] - { - 0, 2, 4 - // todo: vertices of a face must be on lines which share a vertex in the simplex - // that's the pattern, and why it "just works" for the 3-simplex. - // think of it as truncating the embedded simplex - }; + for (var i = 0; i < vlst.Count - 2; i++) + for (var j = i + 1; j < vlst.Count - 1; j++) + for (var k = j + 1; k < vlst.Count; k++) + { + faces.Add(vlst[i]); + faces.Add(vlst[j]); + faces.Add(vlst[k]); - default: - return new int[0]; + edges.Add(vlst[i]); + edges.Add(vlst[j]); + edges.Add(vlst[j]); + edges.Add(vlst[k]); + edges.Add(vlst[k]); + edges.Add(vlst[i]); + } } + + if (faces.Count == 6) + { + for (var i = 0; i < 3; i++) + { + var a = faces[i]; + var b = faces[(1 + i) % 3]; + var c = faces[3 + i]; + var d = faces[3 + (1 + i) % 3]; + + faces.Add(c); + faces.Add(a); + faces.Add(b); + faces.Add(b); + faces.Add(d); + faces.Add(c); + + edges.Add(a); + edges.Add(c); + edges.Add(b); + edges.Add(d); + } + } + + return new Simplex(verts, faces, edges); } } } \ No newline at end of file diff --git a/Tetrahedrons/TetrahedronWindow.cs b/Tetrahedrons/TetrahedronWindow.cs index 124e9da..12521c0 100644 --- a/Tetrahedrons/TetrahedronWindow.cs +++ b/Tetrahedrons/TetrahedronWindow.cs @@ -27,6 +27,7 @@ namespace Tetrahedrons private Simplex _intr; private double _t; + private bool _pause; protected override void OnLoad(EventArgs e) { @@ -62,44 +63,30 @@ namespace Tetrahedrons GL.MatrixMode(MatrixMode.Modelview); GL.LoadMatrix(ref _view); - GL.Begin(PrimitiveType.TriangleStrip); - GL.Color3(1d, 1, 1); - Util.Vertex3(-_a - _b); - Util.Vertex3(-_a + _b); - Util.Vertex3(_a - _b); - Util.Vertex3(_a + _b); + GL.Enable(EnableCap.DepthTest); + + GL.Begin(PrimitiveType.Lines); + GL.Color3(0f, 0, 0); + foreach (var f in _intr.Edjes) + Util.Vertex3(_intr.Verts[f]); GL.End(); - GL.Enable(EnableCap.DepthTest); GL.Begin(PrimitiveType.Triangles); - foreach (var f in _intr.Faces()) + foreach (var f in _intr.Faces) { Util.Color3(_intr.Verts[f]); Util.Vertex3(_intr.Verts[f]); } + GL.End(); GL.Disable(EnableCap.DepthTest); - GL.Begin(PrimitiveType.Points); - foreach (var m in _pent.Verts) - Util.Vertex4(m); - GL.End(); - GL.Begin(PrimitiveType.Lines); - for (var i = 0; i < _pent.Verts.Length; i++) - for (var j = i + 1; j < _pent.Verts.Length; j++) - { - Util.Vertex4(_pent.Verts[i]); - Util.Vertex4(_pent.Verts[j]); - } + foreach (var f in _pent.Edjes) + Util.Vertex4(_pent.Verts[f]); GL.End(); - GL.Begin(PrimitiveType.Points); - foreach (var m in _intr.Verts) - Util.Vertex4(m); - GL.End(); - SwapBuffers(); } @@ -109,22 +96,29 @@ namespace Tetrahedrons _proj3d = Matrix4.CreateOrthographic(6, 6f * Height / Width, -4, 4); + if (_pause) return; + _t += e.Time; - var pln = MVec4D.UnitXy + MVec4D.UnitZw; - + var pln = .5 * MVec4D.UnitXy + MVec4D.UnitZw; var r = MVec4D.Rotor(e.Time / 10, pln.Normalized); for (var i = 0; i < _pent.Verts.Length; i++) - { - var m = _pent.Verts[i]; - m = r * m * !r; - _pent.Verts[i] = m; - } + _pent.Verts[i] |= r; var blade = MVec4D.UnitXyz; -// var pivot = .5 * MVec4D.UnitW; - var pivot = 1.1 * Math.Sin(_t) * MVec4D.UnitW; + var pivot = MVec4D.Zero; + + pivot = .9 * Math.Sin(_t / 5) * MVec4D.UnitW; + _intr = _pent.Intersect(blade, pivot); } + + protected override void OnKeyDown(KeyboardKeyEventArgs e) + { + base.OnKeyDown(e); + + if (e.Key == Key.Space) + _pause = !_pause; + } } } \ No newline at end of file From bd44b4898c262253586aaafaf3c2d1a4d42748d3 Mon Sep 17 00:00:00 2001 From: allem Date: Wed, 4 Apr 2018 18:16:06 -0400 Subject: [PATCH 7/7] bulk rendering --- Tetrahedrons/TetrahedronWindow.cs | 77 +++++++++++++++++++++---------- 1 file changed, 53 insertions(+), 24 deletions(-) diff --git a/Tetrahedrons/TetrahedronWindow.cs b/Tetrahedrons/TetrahedronWindow.cs index 12521c0..a17fdb6 100644 --- a/Tetrahedrons/TetrahedronWindow.cs +++ b/Tetrahedrons/TetrahedronWindow.cs @@ -23,11 +23,12 @@ namespace Tetrahedrons private MVec4D _a = MVec4D.UnitX; private MVec4D _b = MVec4D.UnitY; - private Simplex _pent; - private Simplex _intr; + private Simplex[] _simplexes; + private Simplex[] _intersections; private double _t; private bool _pause; + private bool _wire; protected override void OnLoad(EventArgs e) { @@ -38,12 +39,27 @@ namespace Tetrahedrons _view = Matrix4.LookAt(Vector3.Zero, -new Vector3(.86f, .5f, 1), Vector3.UnitZ); - _pent = new Simplex( // not-quite-regular pentatope - new Vector4d(1, -1, -1, -1), - new Vector4d(-1, 1, -1, -1), - new Vector4d(-1, -1, 1, -1), - new Vector4d(-1, -1, -1, 1), - new Vector4d(1, 1, 1, 1)); + const int n = 5; + _simplexes = new Simplex[n * n * n * n]; + _intersections = new Simplex[_simplexes.Length]; + + var off = (new Vector4d(n - 1, n - 1, n - 1, n - 1)) / 2; + int i = 0; + for (var x = 0; x < n; x++) + for (var y = 0; y < n; y++) + for (var z = 0; z < n; z++) + for (var w = 0; w < n; w++) + { + var d = new Vector4d(x, y, z, w); + _simplexes[i] = new Simplex( // not-quite-regular pentatope + new Vector4d(+.5, -.5, -.5, -.5) + d - off, + new Vector4d(-.5, +.5, -.5, -.5) + d - off, + new Vector4d(-.5, -.5, +.5, -.5) + d - off, + new Vector4d(-.5, -.5, -.5, +.5) + d - off, + new Vector4d(+.5, +.5, +.5, +.5) + d - off); + _intersections[i] = new Simplex(); + i++; + } } protected override void OnRenderFrame(FrameEventArgs e) @@ -67,25 +83,31 @@ namespace Tetrahedrons GL.Begin(PrimitiveType.Lines); GL.Color3(0f, 0, 0); - foreach (var f in _intr.Edjes) - Util.Vertex3(_intr.Verts[f]); + foreach (var intr in _intersections) + foreach (var f in intr.Edjes) + Util.Vertex3(intr.Verts[f]); GL.End(); GL.Begin(PrimitiveType.Triangles); - foreach (var f in _intr.Faces) + foreach (var intr in _intersections) + foreach (var f in intr.Faces) { - Util.Color3(_intr.Verts[f]); - Util.Vertex3(_intr.Verts[f]); + Util.Color3(intr.Verts[f]); + Util.Vertex3(intr.Verts[f]); } GL.End(); GL.Disable(EnableCap.DepthTest); - GL.Begin(PrimitiveType.Lines); - foreach (var f in _pent.Edjes) - Util.Vertex4(_pent.Verts[f]); + if (_wire) + { + GL.Begin(PrimitiveType.Lines); + foreach (var sim in _simplexes) + foreach (var f in sim.Edjes) + Util.Vertex4(sim.Verts[f]); - GL.End(); + GL.End(); + } SwapBuffers(); } @@ -94,23 +116,27 @@ namespace Tetrahedrons { base.OnUpdateFrame(e); - _proj3d = Matrix4.CreateOrthographic(6, 6f * Height / Width, -4, 4); + var w = 10f; + var d = w; + _proj3d = Matrix4.CreateOrthographic(w, w * Height / Width, -d / 2, d / 2); if (_pause) return; _t += e.Time; - var pln = .5 * MVec4D.UnitXy + MVec4D.UnitZw; + var pln = MVec4D.UnitXy + MVec4D.UnitXz + MVec4D.UnitYw; var r = MVec4D.Rotor(e.Time / 10, pln.Normalized); - for (var i = 0; i < _pent.Verts.Length; i++) - _pent.Verts[i] |= r; + foreach (var pent in _simplexes) + for (var i = 0; i < pent.Verts.Length; i++) + pent.Verts[i] |= r; var blade = MVec4D.UnitXyz; var pivot = MVec4D.Zero; - pivot = .9 * Math.Sin(_t / 5) * MVec4D.UnitW; - - _intr = _pent.Intersect(blade, pivot); + for (var i = 0; i < _simplexes.Length; i++) + { + _intersections[i] = _simplexes[i].Intersect(blade, pivot); + } } protected override void OnKeyDown(KeyboardKeyEventArgs e) @@ -119,6 +145,9 @@ namespace Tetrahedrons if (e.Key == Key.Space) _pause = !_pause; + + if (e.Key == Key.Comma) + _wire = !_wire; } } } \ No newline at end of file