Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 3 additions & 16 deletions src/CDT.Core/Predicates.cs
Original file line number Diff line number Diff line change
Expand Up @@ -4,13 +4,6 @@
//
// Geometric predicates — port of Lenthe/Shewchuk adaptive predicates from
// artem-ogre/CDT (predicates.h, predicates::adaptive namespace).
//
// Key design rules matching the C++ template:
// - float inputs: ALL intermediate differences and products computed in float
// first (no premature promotion). Fallback uses decimal exact expansion sign.
// - double inputs: full 3-stage adaptive predicates using Shewchuk/Lenthe
// floating-point expansion arithmetic (all in native double).
// Matches C++ artem-ogre/CDT predicates.h exactly on any IEEE 754 platform.

using System.Runtime.CompilerServices;

Expand All @@ -37,10 +30,7 @@ internal static class Predicates
// Orient2D
// =========================================================================

/// <summary>
/// Adaptive orient2d for <see cref="double"/>. Matches C++
/// <c>predicates::adaptive::orient2d&lt;double&gt;</c> exactly.
/// </summary>
/// <summary>Adaptive orient2d for <see cref="double"/>.</summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static double Orient2D(
double ax, double ay, double bx, double by, double cx, double cy)
Expand Down Expand Up @@ -112,7 +102,7 @@ public static double Orient2D(

/// <summary>
/// Adaptive orient2d for <see cref="float"/>. All fast-path arithmetic in float.
/// Falls back to exact double computation (matching C++).
/// Falls back to exact double computation.
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static float Orient2D(
Expand Down Expand Up @@ -149,10 +139,7 @@ public static float Orient2D(
// InCircle
// =========================================================================

/// <summary>
/// Adaptive incircle for <see cref="double"/>. Matches C++
/// <c>predicates::adaptive::incircle&lt;double&gt;</c> exactly.
/// </summary>
/// <summary>Adaptive incircle for <see cref="double"/>.</summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static double InCircle(
double ax, double ay, double bx, double by,
Expand Down
144 changes: 55 additions & 89 deletions src/CDT.Core/TriangleUtils.cs
Original file line number Diff line number Diff line change
Expand Up @@ -164,50 +164,13 @@ public static int EdgeNeighborFromLocation(PtTriLocation location)
/// to <paramref name="v2"/>.
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static PtLineLocation LocatePointLine(
V2d<double> p, V2d<double> v1, V2d<double> v2, double tolerance = 0.0)
{
double o = Predicates.Orient2D(v1.X, v1.Y, v2.X, v2.Y, p.X, p.Y);
return ClassifyOrientation(o, tolerance);
}

/// <inheritdoc cref="LocatePointLine(V2d{double},V2d{double},V2d{double},double)"/>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static PtLineLocation LocatePointLine(
V2d<float> p, V2d<float> v1, V2d<float> v2, float tolerance = 0f)
{
float o = Predicates.Orient2D(v1.X, v1.Y, v2.X, v2.Y, p.X, p.Y);
return ClassifyOrientation(o, tolerance);
}
public static PtLineLocation LocatePointLine<T>(V2d<T> p, V2d<T> v1, V2d<T> v2, T tolerance = default)
where T : unmanaged, IFloatingPoint<T>
=> ClassifyOrientation(Orient2D(p, v1, v2), tolerance);

/// <summary>Classifies the position of point <paramref name="p"/> within triangle (v1,v2,v3).</summary>
public static PtTriLocation LocatePointTriangle(
V2d<double> p,
V2d<double> v1, V2d<double> v2, V2d<double> v3)
{
PtTriLocation result = PtTriLocation.Inside;

PtLineLocation e = LocatePointLine(p, v1, v2);
if (e == PtLineLocation.Right) return PtTriLocation.Outside;
if (e == PtLineLocation.OnLine) result = PtTriLocation.OnEdge0;

e = LocatePointLine(p, v2, v3);
if (e == PtLineLocation.Right) return PtTriLocation.Outside;
if (e == PtLineLocation.OnLine)
result = result == PtTriLocation.Inside ? PtTriLocation.OnEdge1 : PtTriLocation.OnVertex;

e = LocatePointLine(p, v3, v1);
if (e == PtLineLocation.Right) return PtTriLocation.Outside;
if (e == PtLineLocation.OnLine)
result = result == PtTriLocation.Inside ? PtTriLocation.OnEdge2 : PtTriLocation.OnVertex;

return result;
}

/// <inheritdoc cref="LocatePointTriangle(V2d{double},V2d{double},V2d{double},V2d{double})"/>
public static PtTriLocation LocatePointTriangle(
V2d<float> p,
V2d<float> v1, V2d<float> v2, V2d<float> v3)
public static PtTriLocation LocatePointTriangle<T>(V2d<T> p, V2d<T> v1, V2d<T> v2, V2d<T> v3)
where T : unmanaged, IFloatingPoint<T>
{
PtTriLocation result = PtTriLocation.Inside;

Expand All @@ -232,17 +195,22 @@ public static PtTriLocation LocatePointTriangle(
/// Tests whether point <paramref name="p"/> is inside the circumcircle of triangle (v1,v2,v3).
/// </summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static bool IsInCircumcircle(
V2d<double> p,
V2d<double> v1, V2d<double> v2, V2d<double> v3)
=> Predicates.InCircle(v1.X, v1.Y, v2.X, v2.Y, v3.X, v3.Y, p.X, p.Y) > 0.0;

/// <inheritdoc cref="IsInCircumcircle(V2d{double},V2d{double},V2d{double},V2d{double})"/>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static bool IsInCircumcircle(
V2d<float> p,
V2d<float> v1, V2d<float> v2, V2d<float> v3)
=> Predicates.InCircle(v1.X, v1.Y, v2.X, v2.Y, v3.X, v3.Y, p.X, p.Y) > 0f;
public static bool IsInCircumcircle<T>(V2d<T> p, V2d<T> v1, V2d<T> v2, V2d<T> v3)
where T : unmanaged, IFloatingPoint<T>
{
if (typeof(T) == typeof(double))
return Predicates.InCircle(
Unsafe.As<T, double>(ref Unsafe.AsRef(in v1.X)), Unsafe.As<T, double>(ref Unsafe.AsRef(in v1.Y)),
Unsafe.As<T, double>(ref Unsafe.AsRef(in v2.X)), Unsafe.As<T, double>(ref Unsafe.AsRef(in v2.Y)),
Unsafe.As<T, double>(ref Unsafe.AsRef(in v3.X)), Unsafe.As<T, double>(ref Unsafe.AsRef(in v3.Y)),
Unsafe.As<T, double>(ref Unsafe.AsRef(in p.X)), Unsafe.As<T, double>(ref Unsafe.AsRef(in p.Y))) > 0.0;
else
return Predicates.InCircle(
Unsafe.As<T, float>(ref Unsafe.AsRef(in v1.X)), Unsafe.As<T, float>(ref Unsafe.AsRef(in v1.Y)),
Unsafe.As<T, float>(ref Unsafe.AsRef(in v2.X)), Unsafe.As<T, float>(ref Unsafe.AsRef(in v2.Y)),
Unsafe.As<T, float>(ref Unsafe.AsRef(in v3.X)), Unsafe.As<T, float>(ref Unsafe.AsRef(in v3.Y)),
Unsafe.As<T, float>(ref Unsafe.AsRef(in p.X)), Unsafe.As<T, float>(ref Unsafe.AsRef(in p.Y))) > 0f;
}

/// <summary>Tests whether any pair of triangles (given as vertex triangle lists) share an edge.</summary>
public static bool VerticesShareEdge(List<int> aTris, List<int> bTris)
Expand All @@ -255,45 +223,43 @@ public static bool VerticesShareEdge(List<int> aTris, List<int> bTris)
}

/// <summary>Computes the intersection point of segments (a,b) and (c,d).</summary>
public static V2d<double> IntersectionPosition(
V2d<double> a, V2d<double> b,
V2d<double> c, V2d<double> d)
public static V2d<T> IntersectionPosition<T>(V2d<T> a, V2d<T> b, V2d<T> c, V2d<T> d)
where T : unmanaged, IFloatingPoint<T>
{
double acd = Predicates.Orient2D(c.X, c.Y, d.X, d.Y, a.X, a.Y);
double bcd = Predicates.Orient2D(c.X, c.Y, d.X, d.Y, b.X, b.Y);
double tab = acd / (acd - bcd);

double cab = Predicates.Orient2D(a.X, a.Y, b.X, b.Y, c.X, c.Y);
double dab = Predicates.Orient2D(a.X, a.Y, b.X, b.Y, d.X, d.Y);
double tcd = cab / (cab - dab);

return new V2d<double>(
Math.Abs(a.X - b.X) < Math.Abs(c.X - d.X) ? Lerp(a.X, b.X, tab) : Lerp(c.X, d.X, tcd),
Math.Abs(a.Y - b.Y) < Math.Abs(c.Y - d.Y) ? Lerp(a.Y, b.Y, tab) : Lerp(c.Y, d.Y, tcd));
T acd = Orient2D(a, c, d);
T bcd = Orient2D(b, c, d);
T tab = acd / (acd - bcd);

T cab = Orient2D(c, a, b);
T dab = Orient2D(d, a, b);
T tcd = cab / (cab - dab);

static T Lerp(T x, T y, T t) => (T.One - t) * x + t * y;
return new V2d<T>(
T.Abs(a.X - b.X) < T.Abs(c.X - d.X) ? Lerp(a.X, b.X, tab) : Lerp(c.X, d.X, tcd),
T.Abs(a.Y - b.Y) < T.Abs(c.Y - d.Y) ? Lerp(a.Y, b.Y, tab) : Lerp(c.Y, d.Y, tcd));
}

/// <inheritdoc cref="IntersectionPosition(V2d{double},V2d{double},V2d{double},V2d{double})"/>
public static V2d<float> IntersectionPosition(
V2d<float> a, V2d<float> b,
V2d<float> c, V2d<float> d)
/// <summary>Returns the signed area (orientation) of triangle (v1, v2, p).</summary>
[MethodImpl(MethodImplOptions.AggressiveInlining)]
internal static T Orient2D<T>(V2d<T> p, V2d<T> v1, V2d<T> v2)
where T : unmanaged, IFloatingPoint<T>
{
// Use float-precision Orient2D throughout, matching C++ detail::intersectionPosition<float>.
float acd = Predicates.Orient2D(c.X, c.Y, d.X, d.Y, a.X, a.Y);
float bcd = Predicates.Orient2D(c.X, c.Y, d.X, d.Y, b.X, b.Y);
float tab = acd / (acd - bcd);

float cab = Predicates.Orient2D(a.X, a.Y, b.X, b.Y, c.X, c.Y);
float dab = Predicates.Orient2D(a.X, a.Y, b.X, b.Y, d.X, d.Y);
float tcd = cab / (cab - dab);

return new V2d<float>(
MathF.Abs(a.X - b.X) < MathF.Abs(c.X - d.X) ? LerpF(a.X, b.X, tab) : LerpF(c.X, d.X, tcd),
MathF.Abs(a.Y - b.Y) < MathF.Abs(c.Y - d.Y) ? LerpF(a.Y, b.Y, tab) : LerpF(c.Y, d.Y, tcd));
if (typeof(T) == typeof(double))
{
double r = Predicates.Orient2D(
Unsafe.As<T, double>(ref Unsafe.AsRef(in v1.X)), Unsafe.As<T, double>(ref Unsafe.AsRef(in v1.Y)),
Unsafe.As<T, double>(ref Unsafe.AsRef(in v2.X)), Unsafe.As<T, double>(ref Unsafe.AsRef(in v2.Y)),
Unsafe.As<T, double>(ref Unsafe.AsRef(in p.X)), Unsafe.As<T, double>(ref Unsafe.AsRef(in p.Y)));
return Unsafe.As<double, T>(ref r);
}
else
{
float r = Predicates.Orient2D(
Unsafe.As<T, float>(ref Unsafe.AsRef(in v1.X)), Unsafe.As<T, float>(ref Unsafe.AsRef(in v1.Y)),
Unsafe.As<T, float>(ref Unsafe.AsRef(in v2.X)), Unsafe.As<T, float>(ref Unsafe.AsRef(in v2.Y)),
Unsafe.As<T, float>(ref Unsafe.AsRef(in p.X)), Unsafe.As<T, float>(ref Unsafe.AsRef(in p.Y)));
return Unsafe.As<float, T>(ref r);
}
}

[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static double Lerp(double a, double b, double t) => (1.0 - t) * a + t * b;

[MethodImpl(MethodImplOptions.AggressiveInlining)]
private static float LerpF(float a, float b, float t) => (1f - t) * a + t * b;
}
Loading