CultMath: site cells (power and multiplicative scores, fold, border), Gabriel margin and gate, out-free distance_to_segment

Not yet verified on Yggdrasil; this commit builds only if the next run says so.

Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens 2026-10-16 18:14:50 +02:00
commit a3c3bd9695
6 changed files with 515 additions and 7 deletions

View file

@ -0,0 +1,30 @@
namespace CultMath;
/// <summary>
/// The two nearest sites of a query point under one site score (<see cref="math.power_score"/> or
/// <see cref="math.multiplicative_score"/>), built by folding sites through
/// <see cref="math.site_cells_insert"/>. Each score field follows the
/// <c>float4(gradient.xyz, value.w)</c> layout of a bare value-and-gradient return.
/// </summary>
public struct CultSiteCells
{
/// <summary>(&#8711;F1, F1): the smallest score and its gradient. 1.0e30 in w when no site was inserted.</summary>
public float4 nearest;
/// <summary>(&#8711;F2, F2): the second smallest score and its gradient. 1.0e30 in w when fewer than two sites were inserted.</summary>
public float4 second;
/// <summary>Index of the nearest site as inserted, -1 when none. A float so the struct holds one scalar type, as <see cref="CultCellular.id"/> does; exact below 2^24.</summary>
public float nearestSite;
/// <summary>Index of the second nearest site as inserted, -1 when none.</summary>
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;
}
}

View file

@ -0,0 +1,2 @@
fileFormatVersion: 2
guid: f6479c37f7d349708f54ec7e4d671200

View file

@ -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;
}
/// <summary>
/// The two nearest sites of <paramref name="p"/> by <see cref="power_score"/> (the power diagram), every
/// site visited. A site is <c>float4(position.xyz, weight)</c>; use z = 0 for 2D. Ties keep the lower index first.
/// </summary>
/// <exception cref="ArgumentException">A site weight is not finite.</exception>
public static CultSiteCells power_cells(float3 p, ReadOnlySpan<float4> 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;
}
/// <summary>
/// The two nearest sites of <paramref name="p"/> by <see cref="multiplicative_score"/> (distance divided by
/// weight), every site visited. A site is <c>float4(position.xyz, weight)</c>; use z = 0 for 2D.
/// </summary>
/// <exception cref="ArgumentException">A site weight is not finite and greater than zero.</exception>
public static CultSiteCells multiplicative_cells(float3 p, ReadOnlySpan<float4> 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;
}
/// <summary>
/// Whether sites <paramref name="first"/> and <paramref name="second"/> are neighbours in the Gabriel graph
/// of <paramref name="sites"/> (positions only; weights are ignored): smoothstep(-softness, softness, the
/// least <see cref="gabriel_margin"/> over every other site). 1 when no other site exists.
/// </summary>
/// <exception cref="ArgumentOutOfRangeException">An index is outside the span.</exception>
/// <exception cref="ArgumentException">The indices are equal, or softness is not finite and greater than zero.</exception>
public static float soft_gabriel_gate(ReadOnlySpan<float4> 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);
}
}

View file

@ -0,0 +1,2 @@
fileFormatVersion: 2
guid: 856a5e425b6b45f4a8b3e9a8cd32fe38

View file

@ -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);
}
/// <summary>The distance from <paramref name="point"/> to the segment <paramref name="start"/>-<paramref name="end"/>; a degenerate segment is its start.</summary>
public static float distance_to_segment(float2 point, float2 start, float2 end) =>
distance(point, segment_closest(point, start, end));
/// <summary>The distance from <paramref name="point"/> to the segment <paramref name="start"/>-<paramref name="end"/>; a degenerate segment is its start.</summary>
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)

View file

@ -0,0 +1,314 @@
using Xunit;
using static CultMath.math;
namespace CultMath.Tests;
/// <summary>
/// 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.
/// </summary>
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<float3, float4, float4>, Func<float3, float4, double>)[] { (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<ArgumentException>(() => 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<ArgumentException>(() => 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<ArgumentException>(() => soft_gabriel_gate(three, 1, 1, 0.1f)).ParamName);
Assert.Equal("first", Assert.Throws<ArgumentOutOfRangeException>(() => soft_gabriel_gate(three, -1, 1, 0.1f)).ParamName);
Assert.Equal("second", Assert.Throws<ArgumentOutOfRangeException>(() => 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<ArgumentException>(() => 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));
}
}
}