| | | 1 | | //======================================================================= |
| | | 2 | | // WidePlaneMetrics.cs |
| | | 3 | | //======================================================================= |
| | | 4 | | // MIT License, Copyright (c) 2024-present David Oravsky (mrdav30) |
| | | 5 | | // See LICENSE file in the project root for full license information. |
| | | 6 | | //======================================================================= |
| | | 7 | | |
| | | 8 | | using System; |
| | | 9 | | |
| | | 10 | | namespace FixedMathSharp.Geometry; |
| | | 11 | | |
| | | 12 | | /// <summary>Exact metrics of representable points projected onto an exact rigid plane.</summary> |
| | | 13 | | /// <content>Compact projected metrics and shared exact plane geometry.</content> |
| | | 14 | | internal readonly partial struct WidePlaneMetrics |
| | | 15 | | { |
| | | 16 | | private readonly WideAxis3 normal; |
| | | 17 | | private readonly Signed576 metric; |
| | | 18 | | |
| | | 19 | | /// <summary>Uses a nonzero exact normal with components below 262 bits.</summary> |
| | | 20 | | internal WidePlaneMetrics(in WideAxis3 exactNormal) |
| | | 21 | | { |
| | 3 | 22 | | if (exactNormal.IsZero) throw new ArgumentException("The metric requires a nonzero normal.", nameof(exactNormal) |
| | 1 | 23 | | normal = exactNormal; |
| | 1 | 24 | | metric = normal.SquaredLength; |
| | 1 | 25 | | } |
| | | 26 | | |
| | | 27 | | internal WidePlaneMetrics(FixedTriangle plane, FixedQuaternion rotation) |
| | | 28 | | { |
| | 19 | 29 | | plane.GetExactNormal(out Signed192 x, out Signed192 y, out Signed192 z, out _); |
| | 19 | 30 | | if (x.IsZero && y.IsZero && z.IsZero) |
| | 1 | 31 | | throw new ArgumentException("The metric requires a nondegenerate plane.", nameof(plane)); |
| | 18 | 32 | | var basis = new WideRationalBasis3d(rotation); |
| | 18 | 33 | | if (basis.Denominator.IsZero) |
| | 1 | 34 | | throw new ArgumentException("The metric requires a nonzero rotation.", nameof(rotation)); |
| | 17 | 35 | | normal = WideRigidProjection.TransformLocalAxis(basis, x, y, z); |
| | 17 | 36 | | metric = normal.SquaredLength; |
| | 17 | 37 | | } |
| | | 38 | | |
| | | 39 | | /// <summary>Uses admitted scaled geometry with centered reach at most 4096 units.</summary> |
| | | 40 | | internal WidePlaneMetrics(FixedTriangle plane, in ScaledRigidFrame3d frame) |
| | | 41 | | { |
| | 10 | 42 | | if (frame.IsRigid) { this = new WidePlaneMetrics(plane, frame.Rotation); return; } |
| | 4 | 43 | | WideDirection3d direction = WideTriangleRelations.GetScaledTriangleLocalNormal(plane, frame) |
| | 4 | 44 | | .Transform(frame.GeometryBasis); |
| | 4 | 45 | | normal = new WideAxis3(Signed320.NarrowValue(direction.X), Signed320.NarrowValue(direction.Y), Signed320.NarrowV |
| | 6 | 46 | | if (normal.IsZero) throw new ArgumentException("The metric requires a nondegenerate scaled plane.", nameof(plane |
| | | 47 | | // Centered raw radius <=2^44 gives scaled-edge numerators <=2^109. |
| | | 48 | | // Their cross and one normalized basis give |N|<2^283. Thus M<2^566 |
| | | 49 | | // and full-domain projected area squares <2^830 fit the compact owners. |
| | 2 | 50 | | metric = normal.SquaredLength; |
| | 2 | 51 | | } |
| | | 52 | | |
| | | 53 | | // Raw endpoint differences <2^64, authored normal components <2^130 |
| | | 54 | | // and arbitrary quaternion basis components <2^130 give N<2^262. |
| | | 55 | | // M=N.N<2^526; G(u,v)=M(u.v)-(N.u)(N.v)<2^657; |
| | | 56 | | // T(u,v)=N.(u cross v)<2^393 and T²<2^786. Signed832 retains |
| | | 57 | | // every metric. Ratio comparison uses the existing 26-word products. |
| | | 58 | | |
| | | 59 | | /// <summary>Returns M times squared projected distance in raw coordinate units.</summary> |
| | | 60 | | internal Signed832 GetSpanSquared(Vector3d first, Vector3d second) |
| | | 61 | | { |
| | 70 | 62 | | WideAxis3 edge = Difference(second, first); |
| | 70 | 63 | | return Gram(edge, edge); |
| | | 64 | | } |
| | | 65 | | |
| | | 66 | | /// <summary>Returns 4M times squared projected triangle area in raw coordinate units.</summary> |
| | | 67 | | internal Signed832 GetAreaSquared(Vector3d first, Vector3d second, Vector3d third) |
| | | 68 | | { |
| | 65 | 69 | | Signed576 area = Triple(Difference(second, first), Difference(third, first)); |
| | 65 | 70 | | return WideArithmetic.MultiplySigned576ToSigned832(area, area); |
| | | 71 | | } |
| | | 72 | | |
| | | 73 | | /// <summary> |
| | | 74 | | /// Returns exact squared distance between the projected point and filled |
| | | 75 | | /// projected triangle as a nonnegative numerator over a positive denominator, |
| | | 76 | | /// in raw coordinate units squared. Degenerate projections reduce to edges. |
| | | 77 | | /// </summary> |
| | | 78 | | internal void GetTriangleDistanceSquared(Vector3d point, Vector3d first, Vector3d second, Vector3d third, |
| | | 79 | | out Signed832 numerator, out Signed832 denominator) |
| | | 80 | | { |
| | 606 | 81 | | WideAxis3 edge = Difference(second, first), offset = Difference(point, first); |
| | 303 | 82 | | int orientation = Triple(edge, Difference(third, first)).Sign; |
| | 303 | 83 | | if (orientation != 0 |
| | 303 | 84 | | && orientation * Triple(edge, offset).Sign >= 0 |
| | 303 | 85 | | && orientation * Triple(Difference(third, second), Difference(point, second)).Sign >= 0 |
| | 303 | 86 | | && orientation * Triple(Difference(first, third), Difference(point, third)).Sign >= 0) |
| | | 87 | | { |
| | 160 | 88 | | numerator = default; |
| | 160 | 89 | | denominator = Signed832.ExtendValue(metric); |
| | 160 | 90 | | return; |
| | | 91 | | } |
| | 143 | 92 | | GetSegmentDistanceSquared(point, first, second, out numerator, out denominator); |
| | 143 | 93 | | GetSegmentDistanceSquared(point, second, third, out Signed832 nextNumerator, out Signed832 nextDenominator); |
| | 143 | 94 | | if (CompareRatios(nextNumerator, nextDenominator, numerator, denominator) < 0) |
| | | 95 | | { |
| | 106 | 96 | | numerator = nextNumerator; denominator = nextDenominator; |
| | | 97 | | } |
| | 143 | 98 | | GetSegmentDistanceSquared(point, third, first, out nextNumerator, out nextDenominator); |
| | 143 | 99 | | if (CompareRatios(nextNumerator, nextDenominator, numerator, denominator) < 0) |
| | | 100 | | { |
| | 34 | 101 | | numerator = nextNumerator; denominator = nextDenominator; |
| | | 102 | | } |
| | 143 | 103 | | } |
| | | 104 | | |
| | | 105 | | /// <summary>Compares nonnegative metric ratios with positive denominators.</summary> |
| | | 106 | | internal static int CompareRatios(Signed832 firstNumerator, Signed832 firstDenominator, |
| | | 107 | | Signed832 secondNumerator, Signed832 secondDenominator) => |
| | 436 | 108 | | WideArithmetic.CompareNonNegativeProducts(firstNumerator, secondDenominator, secondNumerator, firstDenominator); |
| | | 109 | | |
| | | 110 | | private void GetSegmentDistanceSquared(Vector3d point, Vector3d first, Vector3d second, |
| | | 111 | | out Signed832 numerator, out Signed832 denominator) |
| | | 112 | | { |
| | 858 | 113 | | WideAxis3 edge = Difference(second, first), offset = Difference(point, first); |
| | 858 | 114 | | Signed832 length = Gram(edge, edge), projection = Gram(offset, edge); |
| | 429 | 115 | | if (length.IsZero || projection.Sign <= 0) |
| | | 116 | | { |
| | 156 | 117 | | numerator = Gram(offset, offset); |
| | 156 | 118 | | denominator = Signed832.ExtendValue(metric); |
| | | 119 | | } |
| | 273 | 120 | | else if (WideArithmetic.SubtractSigned832(projection, length).Sign >= 0) |
| | | 121 | | { |
| | 154 | 122 | | offset = Difference(point, second); |
| | 154 | 123 | | numerator = Gram(offset, offset); |
| | 154 | 124 | | denominator = Signed832.ExtendValue(metric); |
| | | 125 | | } |
| | | 126 | | else |
| | | 127 | | { |
| | | 128 | | // G(v,v)G(e,e)-G(v,e)² = M*T(v,e)². Cancel M before |
| | | 129 | | // forming a numerator, avoiding an unnecessary >832-bit product. |
| | 119 | 130 | | Signed576 area = Triple(offset, edge); |
| | 119 | 131 | | numerator = WideArithmetic.MultiplySigned576ToSigned832(area, area); |
| | 119 | 132 | | denominator = length; |
| | | 133 | | } |
| | 119 | 134 | | } |
| | | 135 | | |
| | | 136 | | private Signed832 Gram(WideAxis3 first, WideAxis3 second) => |
| | 1238 | 137 | | WideArithmetic.SubtractSigned832( |
| | 1238 | 138 | | WideArithmetic.MultiplySigned576ToSigned832(metric, WideAxis3.Dot(first, second)), |
| | 1238 | 139 | | WideArithmetic.MultiplySigned576ToSigned832(WideAxis3.Dot(normal, first), WideAxis3.Dot(normal, second))); |
| | | 140 | | |
| | 1208 | 141 | | private Signed576 Triple(WideAxis3 first, WideAxis3 second) => WideAxis3.Dot(normal, WideAxis3.Cross(first, second)) |
| | | 142 | | |
| | 3916 | 143 | | private static WideAxis3 Difference(Vector3d end, Vector3d start) => new( |
| | 3916 | 144 | | Signed320.ExtendValue(WideArithmetic.SubtractSigned192(Signed192.Raw(end.X), Signed192.Raw(start.X))), |
| | 3916 | 145 | | Signed320.ExtendValue(WideArithmetic.SubtractSigned192(Signed192.Raw(end.Y), Signed192.Raw(start.Y))), |
| | 3916 | 146 | | Signed320.ExtendValue(WideArithmetic.SubtractSigned192(Signed192.Raw(end.Z), Signed192.Raw(start.Z)))); |
| | | 147 | | } |