diff --git a/packages/cultmath/src/CultMath/CultSiteCells.cs b/packages/cultmath/src/CultMath/CultSiteCells.cs new file mode 100644 index 00000000..3e5cea60 --- /dev/null +++ b/packages/cultmath/src/CultMath/CultSiteCells.cs @@ -0,0 +1,30 @@ +namespace CultMath; + +/// +/// The two nearest sites of a query point under one site score ( or +/// ), built by folding sites through +/// . Each score field follows the +/// float4(gradient.xyz, value.w) layout of a bare value-and-gradient return. +/// +public struct CultSiteCells +{ + /// (∇F1, F1): the smallest score and its gradient. 1.0e30 in w when no site was inserted. + public float4 nearest; + + /// (∇F2, F2): the second smallest score and its gradient. 1.0e30 in w when fewer than two sites were inserted. + public float4 second; + + /// Index of the nearest site as inserted, -1 when none. A float so the struct holds one scalar type, as does; exact below 2^24. + public float nearestSite; + + /// Index of the second nearest site as inserted, -1 when none. + public float secondSite; + + public CultSiteCells(float4 nearest, float4 second, float nearestSite, float secondSite) + { + this.nearest = nearest; + this.second = second; + this.nearestSite = nearestSite; + this.secondSite = secondSite; + } +} diff --git a/packages/cultmath/src/CultMath/CultSiteCells.cs.meta b/packages/cultmath/src/CultMath/CultSiteCells.cs.meta new file mode 100644 index 00000000..a3033e23 --- /dev/null +++ b/packages/cultmath/src/CultMath/CultSiteCells.cs.meta @@ -0,0 +1,2 @@ +fileFormatVersion: 2 +guid: f6479c37f7d349708f54ec7e4d671200 diff --git a/packages/cultmath/src/CultMath/math.Sites.cs b/packages/cultmath/src/CultMath/math.Sites.cs new file mode 100644 index 00000000..724415f0 --- /dev/null +++ b/packages/cultmath/src/CultMath/math.Sites.cs @@ -0,0 +1,139 @@ +namespace CultMath; + +public static partial class math +{ + // Sites: the model page's span element float4(position.xyz, weight). A span function is a C# loop over + // a per-site step that has a shader twin (site_cells_empty, then site_cells_insert per site); a shader + // cannot take a span, so shader callers write the same loop over their own buffer. One fold, one rule. + + // Power score of a site: value |p - s|^2 - w^2, gradient 2 (p - s). With this score F2 - F1 is affine in + // the query point, so site_cells_border's first-order distance is the exact distance to the power bisector. + public static float4 power_score(float3 p, float4 site) + { + var offset = p - site.xyz; + return new float4(offset * 2.0f, dot(offset, offset) - site.w * site.w); + } + + // Multiplicative score of a site: value |p - s| / w, gradient (p - s) / (|p - s| w). At p = s the + // gradient is zero, cellular's degenerate-gradient convention. + public static float4 multiplicative_score(float3 p, float4 site) + { + var offset = p - site.xyz; + var r = length(offset); + var gradient = r > 0.0f ? offset / (r * site.w) : new float3(0.0f, 0.0f, 0.0f); + return new float4(gradient, r / site.w); + } + + // The fold's identity: no site yet. 1.0e30 is cellular's sentinel (GLSL ES has no inf literal). + public static CultSiteCells site_cells_empty() => + new(new float4(0.0f, 0.0f, 0.0f, 1.0e30f), new float4(0.0f, 0.0f, 0.0f, 1.0e30f), -1.0f, -1.0f); + + // One fold step. Strict < against F1 (the old F1 moves to second), else strict < against F2, so an exact + // tie keeps the earlier index first. + public static CultSiteCells site_cells_insert(CultSiteCells cells, float4 score, int index) + { + var result = cells; + if (score.w < cells.nearest.w) + { + result.second = cells.nearest; + result.secondSite = cells.nearestSite; + result.nearest = score; + result.nearestSite = (float)index; + } + else if (score.w < cells.second.w) + { + result.second = score; + result.secondSite = (float)index; + } + + return result; + } + + // float4(unit normal, distance) of the border between the two nearest sites: n = grad F2 - grad F1, + // distance (F2 - F1) / |n|, normal n / |n|. First order; exact for power scores. n is the direction in + // which F2 - F1 grows, into the nearest site's cell. (0, 0, 0, 1.0e30) when there is no second site or + // |n| = 0. + public static float4 site_cells_border(CultSiteCells cells) + { + var n = cells.second.xyz - cells.nearest.xyz; + var len = length(n); + if (cells.secondSite < 0.0f || len <= 0.0f) + return new float4(0.0f, 0.0f, 0.0f, 1.0e30f); + return new float4(n / len, (cells.second.w - cells.nearest.w) / len); + } + + // Gabriel margin of the pair (a, b) against a third site: (|other - mid| - r) / r with mid = (a + b) / 2 + // and r = |a - b| / 2. 0 on the diameter circle, negative inside it, -1 when a and b coincide. + public static float gabriel_margin(float3 a, float3 b, float3 other) + { + var r = distance(a, b) * 0.5f; + if (r <= 0.0f) + return -1.0f; + return (distance(other, (a + b) * 0.5f) - r) / r; + } + + /// + /// The two nearest sites of by (the power diagram), every + /// site visited. A site is float4(position.xyz, weight); use z = 0 for 2D. Ties keep the lower index first. + /// + /// A site weight is not finite. + public static CultSiteCells power_cells(float3 p, ReadOnlySpan sites) + { + var cells = site_cells_empty(); + for (var i = 0; i < sites.Length; i++) + { + if (!float.IsFinite(sites[i].w)) + throw new ArgumentException("Site weights must be finite.", nameof(sites)); + cells = site_cells_insert(cells, power_score(p, sites[i]), i); + } + + return cells; + } + + /// + /// The two nearest sites of by (distance divided by + /// weight), every site visited. A site is float4(position.xyz, weight); use z = 0 for 2D. + /// + /// A site weight is not finite and greater than zero. + public static CultSiteCells multiplicative_cells(float3 p, ReadOnlySpan sites) + { + var cells = site_cells_empty(); + for (var i = 0; i < sites.Length; i++) + { + if (!(float.IsFinite(sites[i].w) && sites[i].w > 0.0f)) + throw new ArgumentException("Site weights must be finite and greater than zero.", nameof(sites)); + cells = site_cells_insert(cells, multiplicative_score(p, sites[i]), i); + } + + return cells; + } + + /// + /// Whether sites and are neighbours in the Gabriel graph + /// of (positions only; weights are ignored): smoothstep(-softness, softness, the + /// least over every other site). 1 when no other site exists. + /// + /// An index is outside the span. + /// The indices are equal, or softness is not finite and greater than zero. + public static float soft_gabriel_gate(ReadOnlySpan sites, int first, int second, float softness) + { + if (first < 0 || first >= sites.Length) + throw new ArgumentOutOfRangeException(nameof(first), "Index is outside the site span."); + if (second < 0 || second >= sites.Length) + throw new ArgumentOutOfRangeException(nameof(second), "Index is outside the site span."); + if (first == second) + throw new ArgumentException("The two indices must differ.", nameof(second)); + if (!(float.IsFinite(softness) && softness > 0.0f)) + throw new ArgumentException("Softness must be finite and greater than zero.", nameof(softness)); + + var margin = float.PositiveInfinity; + for (var k = 0; k < sites.Length; k++) + { + if (k == first || k == second) + continue; + margin = min(margin, gabriel_margin(sites[first].xyz, sites[second].xyz, sites[k].xyz)); + } + + return smoothstep(-softness, softness, margin); + } +} diff --git a/packages/cultmath/src/CultMath/math.Sites.cs.meta b/packages/cultmath/src/CultMath/math.Sites.cs.meta new file mode 100644 index 00000000..f9956189 --- /dev/null +++ b/packages/cultmath/src/CultMath/math.Sites.cs.meta @@ -0,0 +1,2 @@ +fileFormatVersion: 2 +guid: 856a5e425b6b45f4a8b3e9a8cd32fe38 diff --git a/packages/cultmath/src/CultMath/math.cs b/packages/cultmath/src/CultMath/math.cs index a5a30a09..59197262 100644 --- a/packages/cultmath/src/CultMath/math.cs +++ b/packages/cultmath/src/CultMath/math.cs @@ -317,18 +317,39 @@ public static partial class math } public static float distance_to_segment(float2 point, float2 start, float2 end, out float2 closest) + { + closest = segment_closest(point, start, end); + return distance(point, closest); + } + + /// The distance from to the segment -; a degenerate segment is its start. + public static float distance_to_segment(float2 point, float2 start, float2 end) => + distance(point, segment_closest(point, start, end)); + + /// The distance from to the segment -; a degenerate segment is its start. + public static float distance_to_segment(float3 point, float3 start, float3 end) => + distance(point, segment_closest(point, start, end)); + + private static float2 segment_closest(float2 point_value, float2 start, float2 end) { var segment = end - start; var segmentLengthSq = lengthsq(segment); if (segmentLengthSq <= 0.0f) - { - closest = start; - return distance(point, start); - } + return start; - var t = saturate(dot(point - start, segment) / segmentLengthSq); - closest = start + segment * t; - return distance(point, closest); + var t = saturate(dot(point_value - start, segment) / segmentLengthSq); + return start + segment * t; + } + + private static float3 segment_closest(float3 point_value, float3 start, float3 end) + { + var segment = end - start; + var segmentLengthSq = lengthsq(segment); + if (segmentLengthSq <= 0.0f) + return start; + + var t = saturate(dot(point_value - start, segment) / segmentLengthSq); + return start + segment * t; } public static float catmullrom(float p0, float p1, float p2, float p3, float t) diff --git a/packages/cultmath/tests/CultMath.Tests/SiteCellsTests.cs b/packages/cultmath/tests/CultMath.Tests/SiteCellsTests.cs new file mode 100644 index 00000000..12859948 --- /dev/null +++ b/packages/cultmath/tests/CultMath.Tests/SiteCellsTests.cs @@ -0,0 +1,314 @@ +using Xunit; +using static CultMath.math; + +namespace CultMath.Tests; + +/// +/// The site cells (power and multiplicative scores, the fold, the border), the Gabriel margin and gate, and the +/// out-free segment distance. Every reference is written here in double from the site data, never by calling the +/// function under test. +/// +public sealed class SiteCellsTests +{ + private static float4[] RandomSites(System.Random r, int max) + { + var count = r.Next(max + 1); + var sites = new float4[count]; + for (var i = 0; i < count; i++) + { + sites[i] = i > 0 && r.Next(4) == 0 + ? sites[r.Next(i)] + : new float4(Span(r, 10.0f), Span(r, 10.0f), Span(r, 10.0f), 0.5f + 2.5f * r.NextSingle()); + } + + return sites; + } + + private static float Span(System.Random r, float extent) => (r.NextSingle() * 2.0f - 1.0f) * extent; + + private static float3 Point(System.Random r) => new(Span(r, 10.0f), Span(r, 10.0f), Span(r, 10.0f)); + + private static double PowerReference(float3 p, float4 s) + { + double dx = (double)p.x - s.x, dy = (double)p.y - s.y, dz = (double)p.z - s.z; + return dx * dx + dy * dy + dz * dz - (double)s.w * s.w; + } + + private static double MultiplicativeReference(float3 p, float4 s) + { + double dx = (double)p.x - s.x, dy = (double)p.y - s.y, dz = (double)p.z - s.z; + return Math.Sqrt(dx * dx + dy * dy + dz * dz) / s.w; + } + + // Checks cells against the stable (score, index) order of a double reference. A pair whose reference scores + // differ by less than a float's resolution may legitimately order either way, so such a pair is skipped; an + // exact tie (a duplicated site) is compared. + private static int CheckFold(CultSiteCells cells, double[] reference) + { + var order = Enumerable.Range(0, reference.Length).OrderBy(i => reference[i]).ThenBy(i => i).ToArray(); + if (order.Length < 2) + { + Assert.Equal(-1.0f, cells.secondSite); + Assert.Equal(1.0e30f, cells.second.w); + } + + if (order.Length < 1) + { + Assert.Equal(-1.0f, cells.nearestSite); + Assert.Equal(1.0e30f, cells.nearest.w); + return 0; + } + + bool Clear(int a, int b) => reference[a] == reference[b] || Math.Abs(reference[b] - reference[a]) > 1.0e-4 * (1.0 + Math.Abs(reference[b])); + if (order.Length >= 2 && !Clear(order[0], order[1])) + return 0; + Assert.Equal((float)order[0], cells.nearestSite); + if (order.Length < 2 || (order.Length > 2 && !Clear(order[1], order[2]))) + return 1; + Assert.Equal((float)order[1], cells.secondSite); + return 1; + } + + [Fact] + public void PowerAndMultiplicativeCellsAreTheTwoSmallestScoresInIndexOrder() + { + var r = new System.Random(0x517E); + var checkedSpans = 0; + for (var n = 0; n < 4000; n++) + { + var sites = RandomSites(r, 12); + var p = Point(r); + checkedSpans += CheckFold(power_cells(p, sites), sites.Select(s => PowerReference(p, s)).ToArray()); + checkedSpans += CheckFold(multiplicative_cells(p, sites), sites.Select(s => MultiplicativeReference(p, s)).ToArray()); + } + + Assert.True(checkedSpans > 5000, $"only {checkedSpans} spans were comparable"); + } + + [Fact] + public void ExactTiesKeepTheLowerIndexFirst() + { + var site = new float4(1.0f, 2.0f, 3.0f, 1.5f); + var sites = new[] { site, site, site }; + var cells = power_cells(new float3(0.0f, 0.0f, 0.0f), sites); + Assert.Equal(0.0f, cells.nearestSite); + Assert.Equal(1.0f, cells.secondSite); + } + + // TN2: F2 - F1 is affine for power scores, so the border distance is the distance to the bisector plane + // {x : |x - s1|^2 - w1^2 = |x - s2|^2 - w2^2}, computed here in double from the two winning sites. + [Fact] + public void PowerBorderIsTheDistanceToThePowerBisector() + { + var r = new System.Random(0xB0DE); + var compared = 0; + for (var n = 0; n < 3000; n++) + { + var sites = RandomSites(r, 8); + if (sites.Length < 3) + continue; + var p = Point(r); + var reference = sites.Select(s => PowerReference(p, s)).ToArray(); + var order = Enumerable.Range(0, sites.Length).OrderBy(i => reference[i]).ThenBy(i => i).ToArray(); + if (reference[order[1]] - reference[order[0]] < 1.0e-2 || reference[order[2]] - reference[order[1]] < 1.0e-2) + continue; + var (a, b) = (sites[order[0]], sites[order[1]]); + double nx = 2.0 * ((double)a.x - b.x), ny = 2.0 * ((double)a.y - b.y), nz = 2.0 * ((double)a.z - b.z); + var norm = Math.Sqrt(nx * nx + ny * ny + nz * nz); + if (norm < 1.0e-2) + continue; + var expected = (reference[order[1]] - reference[order[0]]) / norm; + + var border = site_cells_border(power_cells(p, sites)); + Assert.InRange(border.w, expected - 1.0e-4 * expected - 1.0e-3, expected + 1.0e-4 * expected + 1.0e-3); + // F2 - F1 grows toward the nearest site: the normal is (s1 - s2) / |s1 - s2|. + Assert.InRange(border.x, nx / norm - 1.0e-3, nx / norm + 1.0e-3); + Assert.InRange(border.y, ny / norm - 1.0e-3, ny / norm + 1.0e-3); + Assert.InRange(border.z, nz / norm - 1.0e-3, nz / norm + 1.0e-3); + compared++; + } + + Assert.True(compared > 1000, $"only {compared} borders were comparable"); + } + + [Fact] + public void BorderIsEmptyWithoutASecondSite() + { + var one = power_cells(new float3(1.0f, 2.0f, 3.0f), new[] { new float4(0.0f, 0.0f, 0.0f, 1.0f) }); + Assert.Equal(new float4(0.0f, 0.0f, 0.0f, 1.0e30f), site_cells_border(one)); + Assert.Equal(new float4(0.0f, 0.0f, 0.0f, 1.0e30f), site_cells_border(site_cells_empty())); + } + + // Two sites with weights w1 != w2: the border |p - s1| / w1 = |p - s2| / w2 is the Apollonius circle, centre + // (w1^2 s2 - w2^2 s1) / (w1^2 - w2^2), radius w1 w2 |s1 - s2| / |w1^2 - w2^2| (built in double). On it the + // multiplicative border distance vanishes; off it the distance is nonzero. + [Fact] + public void MultiplicativeBorderVanishesOnTheApolloniusCircle() + { + var r = new System.Random(0xA9011); + for (var n = 0; n < 300; n++) + { + var s1 = new float4(Span(r, 5.0f), Span(r, 5.0f), 0.0f, 1.0f + 2.0f * r.NextSingle()); + var s2 = new float4(Span(r, 5.0f), Span(r, 5.0f), 0.0f, 1.0f + 2.0f * r.NextSingle()); + double w1 = s1.w, w2 = s2.w, k = w1 * w1 - w2 * w2; + double sep = Math.Sqrt(Math.Pow((double)s1.x - s2.x, 2) + Math.Pow((double)s1.y - s2.y, 2)); + if (Math.Abs(k) < 0.5 || sep < 1.0) + continue; + double cx = (w1 * w1 * s2.x - w2 * w2 * s1.x) / k, cy = (w1 * w1 * s2.y - w2 * w2 * s1.y) / k; + var radius = w1 * w2 * sep / Math.Abs(k); + if (radius > 50.0) + continue; + + var sites = new[] { s1, s2 }; + var angle = r.NextDouble() * 2.0 * Math.PI; + var on = new float3((float)(cx + radius * Math.Cos(angle)), (float)(cy + radius * Math.Sin(angle)), 0.0f); + var onBorder = site_cells_border(multiplicative_cells(on, sites)); + Assert.True(Math.Abs(onBorder.w) < 1.0e-3 * radius + 1.0e-3, $"on the circle: {onBorder.w}"); + + var off = new float3((float)(cx + 1.5 * radius * Math.Cos(angle)), (float)(cy + 1.5 * radius * Math.Sin(angle)), 0.0f); + var offBorder = site_cells_border(multiplicative_cells(off, sites)); + Assert.True(Math.Abs(offBorder.w) > 1.0e-2 * radius, $"off the circle: {offBorder.w}"); + } + } + + // Central differences in double of each score's value against the returned gradient. + [Fact] + public void ScoreGradientsAreAnalytic() + { + var r = new System.Random(0x6AAD); + const double h = 1.0e-3; + for (var n = 0; n < 500; n++) + { + var site = RandomSites(r, 1).FirstOrDefault(new float4(1.0f, -2.0f, 0.5f, 1.5f)); + var p = Point(r); + if (Math.Abs(p.x - site.x) + Math.Abs(p.y - site.y) + Math.Abs(p.z - site.z) < 1.0f) + continue; + foreach (var (score, reference) in new (Func, Func)[] { (power_score, PowerReference), (multiplicative_score, MultiplicativeReference) }) + { + var got = score(p, site); + var axes = new[] { new float3((float)h, 0.0f, 0.0f), new float3(0.0f, (float)h, 0.0f), new float3(0.0f, 0.0f, (float)h) }; + var expected = axes.Select(a => (reference(p + a, site) - reference(p - a, site)) / (2.0 * h)).ToArray(); + var scale = 1.0 + Math.Abs(expected[0]) + Math.Abs(expected[1]) + Math.Abs(expected[2]); + Assert.True(Math.Abs(got.x - expected[0]) < 2.0e-3 * scale, "x"); + Assert.True(Math.Abs(got.y - expected[1]) < 2.0e-3 * scale, "y"); + Assert.True(Math.Abs(got.z - expected[2]) < 2.0e-3 * scale, "z"); + Assert.Equal(reference(p, site), got.w, 1.0e-3 * (1.0 + Math.Abs(got.w))); + } + } + } + + [Fact] + public void MultiplicativeGradientIsZeroAtTheSite() + { + var site = new float4(1.0f, 2.0f, 3.0f, 2.0f); + var score = multiplicative_score(new float3(1.0f, 2.0f, 3.0f), site); + Assert.Equal(new float4(0.0f, 0.0f, 0.0f, 0.0f), score); + } + + private static double MarginReference(float4 a, float4 b, float4 o) + { + double mx = ((double)a.x + b.x) / 2, my = ((double)a.y + b.y) / 2, mz = ((double)a.z + b.z) / 2; + var rr = Math.Sqrt(Math.Pow((double)a.x - b.x, 2) + Math.Pow((double)a.y - b.y, 2) + Math.Pow((double)a.z - b.z, 2)) / 2; + var d = Math.Sqrt(Math.Pow(o.x - mx, 2) + Math.Pow(o.y - my, 2) + Math.Pow(o.z - mz, 2)); + return (d - rr) / rr; + } + + private static double SmoothstepReference(double lo, double hi, double x) + { + var t = Math.Clamp((x - lo) / (hi - lo), 0.0, 1.0); + return t * t * (3.0 - 2.0 * t); + } + + [Fact] + public void SoftGabrielGateFollowsTheEmptyDiscMargin() + { + var a = new float4(-2.0f, 0.0f, 0.0f, 1.0f); + var b = new float4(2.0f, 0.0f, 0.0f, 1.0f); + // On the diameter circle (radius 2 about the origin): margin 0, gate 0.5. At half radius: margin -0.5, gate 0. + Assert.Equal(0.5f, soft_gabriel_gate(new[] { a, b, new float4(0.0f, 2.0f, 0.0f, 1.0f) }, 0, 1, 0.15f), 6); + Assert.Equal(0.0f, soft_gabriel_gate(new[] { a, b, new float4(0.0f, 1.0f, 0.0f, 1.0f) }, 0, 1, 0.15f)); + Assert.Equal(1.0f, soft_gabriel_gate(new[] { a, b }, 0, 1, 0.15f)); + + var r = new System.Random(0x6ABE); + for (var n = 0; n < 500; n++) + { + var sites = RandomSites(r, 9); + if (sites.Length < 3 || sites[0].Equals(sites[1])) + continue; + var first = r.Next(sites.Length); + var second = (first + 1 + r.Next(sites.Length - 1)) % sites.Length; + var softness = 0.05f + 0.5f * r.NextSingle(); + var margin = Enumerable.Range(0, sites.Length).Where(k => k != first && k != second).Min(k => MarginReference(sites[first], sites[second], sites[k])); + if (!double.IsFinite(margin)) + continue; + Assert.Equal(SmoothstepReference(-softness, softness, margin), soft_gabriel_gate(sites, first, second, softness), 1.0e-3); + } + + // Continuity: a tiny radial step of a third site moves the gate by a small share of the softness band. + var prev = soft_gabriel_gate(new[] { a, b, new float4(0.0f, 1.8f, 0.0f, 1.0f) }, 0, 1, 0.3f); + for (var y = 1.8f; y < 2.2f; y += 0.001f) + { + var gate = soft_gabriel_gate(new[] { a, b, new float4(0.0f, y, 0.0f, 1.0f) }, 0, 1, 0.3f); + Assert.True(Math.Abs(gate - prev) < 0.01f, $"step at y={y}"); + prev = gate; + } + } + + [Fact] + public void SpanFunctionsRefuseBadInputsWithoutEchoingThem() + { + var good = new float4(0.0f, 0.0f, 0.0f, 1.0f); + var marker = 12345.678f; + foreach (var weight in new[] { 0.0f, -marker, float.NaN, float.PositiveInfinity }) + { + var e = Assert.Throws(() => multiplicative_cells(default, new[] { good, new float4(1.0f, 1.0f, 1.0f, weight) })); + Assert.Equal("sites", e.ParamName); + Assert.DoesNotContain("12345", e.Message); + } + + foreach (var weight in new[] { float.NaN, float.PositiveInfinity, float.NegativeInfinity }) + Assert.Equal("sites", Assert.Throws(() => power_cells(default, new[] { good, new float4(1.0f, 1.0f, 1.0f, weight) })).ParamName); + + var three = new[] { good, new float4(1.0f, 0.0f, 0.0f, 1.0f), new float4(0.0f, 1.0f, 0.0f, 1.0f) }; + Assert.Equal("second", Assert.Throws(() => soft_gabriel_gate(three, 1, 1, 0.1f)).ParamName); + Assert.Equal("first", Assert.Throws(() => soft_gabriel_gate(three, -1, 1, 0.1f)).ParamName); + Assert.Equal("second", Assert.Throws(() => soft_gabriel_gate(three, 0, 3, 0.1f)).ParamName); + foreach (var softness in new[] { 0.0f, -marker, float.NaN, float.PositiveInfinity }) + { + var e = Assert.Throws(() => soft_gabriel_gate(three, 0, 1, softness)); + Assert.Equal("softness", e.ParamName); + Assert.DoesNotContain("12345", e.Message); + } + } + + [Fact] + public void OutFreeDistanceToSegmentEqualsTheOutForm() + { + var r = new System.Random(0x5E61); + for (var n = 0; n < 2000; n++) + { + var p = new float2(Span(r, 5.0f), Span(r, 5.0f)); + var start = new float2(Span(r, 5.0f), Span(r, 5.0f)); + var end = r.Next(8) == 0 ? start : new float2(Span(r, 5.0f), Span(r, 5.0f)); + var outForm = distance_to_segment(p, start, end, out var closest); + Assert.Equal(BitConverter.SingleToInt32Bits(outForm), BitConverter.SingleToInt32Bits(distance_to_segment(p, start, end))); + Assert.Equal(distance(p, closest), outForm); + + var p3 = new float3(Span(r, 5.0f), Span(r, 5.0f), Span(r, 5.0f)); + var a3 = new float3(Span(r, 5.0f), Span(r, 5.0f), Span(r, 5.0f)); + var b3 = r.Next(8) == 0 ? a3 : new float3(Span(r, 5.0f), Span(r, 5.0f), Span(r, 5.0f)); + var got = distance_to_segment(p3, a3, b3); + if (a3.Equals(b3)) + { + Assert.Equal(distance(p3, a3), got); + continue; + } + + double ex = (double)b3.x - a3.x, ey = (double)b3.y - a3.y, ez = (double)b3.z - a3.z; + double px = (double)p3.x - a3.x, py = (double)p3.y - a3.y, pz = (double)p3.z - a3.z; + var t = Math.Clamp((px * ex + py * ey + pz * ez) / (ex * ex + ey * ey + ez * ez), 0.0, 1.0); + var expected = Math.Sqrt(Math.Pow(px - t * ex, 2) + Math.Pow(py - t * ey, 2) + Math.Pow(pz - t * ez, 2)); + Assert.Equal(expected, got, 1.0e-4 * (1.0 + expected)); + } + } +}