Port upstream CDT refinement and topology updates

LICENSE-predicates.txt

@@ -0,0 +1,27 @@
+Copyright (c) 2019, William C. Lenthe
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+   contributors may be used to endorse or promote products derived from
+   this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

README.md

@@ -6,6 +6,7 @@
 A C# port of the [artem-ogre/CDT](https://github.com/artem-ogre/CDT) Constrained Delaunay Triangulation library.
 
 > **Credits:** This library is a C# port of the excellent [CDT C++ library](https://github.com/artem-ogre/CDT) by Artem Amirkhanov and contributors, licensed under MPL 2.0.
+> The adaptive geometric predicates are ported from `predicates.h` by William C. Lenthe under BSD 3-Clause; see [`LICENSE-predicates.txt`](LICENSE-predicates.txt).
 > For full algorithm documentation, research references, and in-depth API documentation please refer to the [original C++ repository](https://github.com/artem-ogre/CDT) and its [online documentation](https://artem-ogre.github.io/CDT/).
 
 ## What is CDT?
@@ -16,6 +17,7 @@ CDT is a library for generating **Constrained** and **Conforming** Delaunay Tria
 
 - **Constrained Delaunay Triangulation** — force specific edges into the triangulation
 - **Conforming Delaunay Triangulation** — split edges and add Steiner points until constraint edges are present in the triangulation
+- **Delaunay refinement** — improve triangle quality using minimum-angle and maximum-area criteria
 - **Convex-hull** triangulation — triangulate all points without any constraints
 - **Automatic hole detection** — use `EraseOuterTrianglesAndHoles` to remove outer regions and holes based on even–odd winding depth
 - **Robust geometric predicates** — numerically stable orientation and in-circle tests

src/CDT.Core/CDT.Core.csproj

@@ -23,7 +23,7 @@
     <PackageProjectUrl>https://github.com/MichaCo/CDT.NET</PackageProjectUrl>
     <RepositoryUrl>https://github.com/MichaCo/CDT.NET</RepositoryUrl>
     <RepositoryType>git</RepositoryType>
-    <PackageLicenseExpression>MPL-2.0</PackageLicenseExpression>
+    <PackageLicenseExpression>MPL-2.0 AND BSD-3-Clause</PackageLicenseExpression>
     <PackageReadmeFile>README.md</PackageReadmeFile>
     <PublishRepositoryUrl>true</PublishRepositoryUrl>
     <EmbedUntrackedSources>true</EmbedUntrackedSources>
@@ -33,7 +33,9 @@
   </PropertyGroup>
 
   <ItemGroup>
+    <None Include="..\..\LICENSE" Pack="true" PackagePath="\" />
     <None Include="..\..\README.md" Pack="true" PackagePath="\" />
+    <None Include="..\..\LICENSE-predicates.txt" Pack="true" PackagePath="\" />
   </ItemGroup>
 
   <ItemGroup>

src/CDT.Core/CdtUtils.cs

@@ -541,6 +541,70 @@ public static bool IsInCircumcircle<T>(V2d<T> p, V2d<T> v1, V2d<T> v2, V2d<T> v3
                 Unsafe.As<T, float>(ref Unsafe.AsRef(in p.X)), Unsafe.As<T, float>(ref Unsafe.AsRef(in p.Y))) > 0f;
     }
 
+    /// <summary>
+    /// Returns whether <paramref name="p"/> lies strictly inside the circle
+    /// whose diameter is the segment from <paramref name="a"/> to <paramref name="b"/>.
+    /// </summary>
+    public static bool IsEncroachingOnEdge<T>(V2d<T> p, V2d<T> a, V2d<T> b)
+        where T : unmanaged, IFloatingPoint<T>
+        => PredicatesAdaptive.InDiamCircle(
+            double.CreateChecked(a.X), double.CreateChecked(a.Y),
+            double.CreateChecked(b.X), double.CreateChecked(b.Y),
+            double.CreateChecked(p.X), double.CreateChecked(p.Y)) > 0.0;
+
+    /// <summary>Computes the area of triangle <paramref name="a"/>, <paramref name="b"/>, <paramref name="c"/>.</summary>
+    public static T Area<T>(V2d<T> a, V2d<T> b, V2d<T> c)
+        where T : unmanaged, IFloatingPoint<T>
+        => T.Abs(Orient2D(a, b, c)) / (T.One + T.One);
+
+    /// <summary>Computes the circumcenter of a non-degenerate triangle.</summary>
+    public static V2d<T> Circumcenter<T>(V2d<T> a, V2d<T> b, V2d<T> c)
+        where T : unmanaged, IFloatingPoint<T>
+    {
+        T two = T.One + T.One;
+        T denom = two * Orient2D(a, b, c);
+        if (denom == T.Zero) throw new ArgumentException("Cannot compute the circumcenter of a degenerate triangle.");
+
+        T aLenSq = DistanceSquared(a, c);
+        T bLenSq = DistanceSquared(b, c);
+        a = new V2d<T>(a.X - c.X, a.Y - c.Y);
+        b = new V2d<T>(b.X - c.X, b.Y - c.Y);
+        return new V2d<T>(
+            c.X + (b.Y * aLenSq - a.Y * bLenSq) / denom,
+            c.Y + (a.X * bLenSq - b.X * aLenSq) / denom);
+    }
+
+    /// <summary>Computes the smallest interior angle of a non-degenerate triangle, in radians.</summary>
+    public static T SmallestAngle<T>(V2d<T> a, V2d<T> b, V2d<T> c)
+        where T : unmanaged, IFloatingPoint<T>, IRootFunctions<T>
+    {
+        T sideA = Distance(a, b);
+        T sideB = Distance(b, c);
+        if (sideA > sideB) (sideA, sideB) = (sideB, sideA);
+        sideA = T.Max(sideA, Distance(a, c));
+        T sine = (T.One + T.One) * Area(a, b, c) / sideA / sideB;
+        return T.CreateChecked(Math.Asin(double.CreateChecked(sine)));
+    }
+
+    /// <summary>Converts an angle in degrees to radians.</summary>
+    public static T DegToRad<T>(T degrees)
+        where T : IFloatingPoint<T>
+        => degrees / T.CreateChecked(180) * T.Pi;
+
+    internal static bool IsAngleAtApexSmall<T>(V2d<T> apex, V2d<T> a, V2d<T> b)
+        where T : unmanaged, IFloatingPoint<T>, IRootFunctions<T>
+    {
+        T ux = a.X - apex.X, uy = a.Y - apex.Y;
+        T wx = b.X - apex.X, wy = b.Y - apex.Y;
+        T dot = ux * wx + uy * wy;
+        return dot > T.Zero &&
+               (T.One + T.One) * dot > Distance(apex, a) * Distance(apex, b);
+    }
+
+    internal static T Lerp<T>(T a, T b, T t)
+        where T : IFloatingPoint<T>
+        => (T.One - t) * a + t * b;
+
     /// <summary>
     /// Computes the squared Euclidean distance between two points.
     /// </summary>

src/CDT.Core/Predicates/PredicatesAdaptive.cs

@@ -4,6 +4,7 @@
 //
 // Geometric predicates — port of Lenthe/Shewchuk adaptive predicates from
 // artem-ogre/CDT (predicates.h, predicates::adaptive namespace).
+// Copyright (c) 2019, William C. Lenthe. See LICENSE-predicates.txt.
 
 using System.Runtime.CompilerServices;
 using System.Runtime.InteropServices;
@@ -334,6 +335,58 @@ public static float InCircle(
         return PredicatesExact.InCircle(ax, ay, bx, by, cx, cy, dx, dy);
     }
 
+    /// <summary>
+    /// Returns positive when point c lies strictly inside the circle with
+    /// diameter ab, zero when it lies on the circle, and negative when outside.
+    /// </summary>
+    public static double InDiamCircle(
+        double ax, double ay, double bx, double by, double cx, double cy)
+    {
+        double acx = ax - cx, acy = ay - cy;
+        double bcx = bx - cx, bcy = by - cy;
+        double nbcx = -bcx;
+        double detleft = acx * nbcx;
+        double detright = acy * bcy;
+        double det = detleft - detright;
+
+        if ((detleft < 0.0) != (detright < 0.0)) return det;
+        if (detleft == 0.0 || detright == 0.0) return det;
+
+        double detsum = Math.Abs(detleft + detright);
+        if (Math.Abs(det) >= CcwBoundAD * detsum) return det;
+
+        Span<double> b = stackalloc double[4];
+        int bLen = TwoTwoDiff(acx, nbcx, acy, bcy, b);
+        det = Estimate(b, bLen);
+        if (Math.Abs(det) >= CcwBoundBD * detsum) return det;
+
+        double acxtail = MinusTail(ax, cx, acx);
+        double acytail = MinusTail(ay, cy, acy);
+        double bcxtail = MinusTail(bx, cx, bcx);
+        double bcytail = MinusTail(by, cy, bcy);
+        double nbcxtail = -bcxtail;
+        if (acxtail == 0.0 && nbcxtail == 0.0 && acytail == 0.0 && bcytail == 0.0)
+            return det;
+
+        double errbound = CcwBoundCD * detsum + ResultErrBound * Math.Abs(det);
+        det += (acx * nbcxtail + nbcx * acxtail) - (acy * bcytail + bcy * acytail);
+        if (Math.Abs(det) >= errbound) return det;
+
+        Span<double> s1 = stackalloc double[4];
+        int s1Len = TwoTwoDiff(acxtail, nbcx, acytail, bcy, s1);
+        Span<double> s2 = stackalloc double[4];
+        int s2Len = TwoTwoDiff(acx, nbcxtail, acy, bcytail, s2);
+        Span<double> s3 = stackalloc double[4];
+        int s3Len = TwoTwoDiff(acxtail, nbcxtail, acytail, bcytail, s3);
+        Span<double> t1 = stackalloc double[8];
+        int t1Len = ExpansionSum(b, bLen, s1, s1Len, t1);
+        Span<double> t2 = stackalloc double[12];
+        int t2Len = ExpansionSum(t1, t1Len, s2, s2Len, t2);
+        Span<double> d = stackalloc double[16];
+        int dLen = ExpansionSum(t2, t2Len, s3, s3Len, d);
+        return MostSignificant(d, dLen);
+    }
+
     // =========================================================================
     // Shewchuk/Lenthe floating-point expansion primitives (internal)
     // =========================================================================

src/CDT.Core/TopologyVerifier.cs

@@ -53,5 +53,29 @@ public static bool VerifyTopology<T>(Triangulation<T> cdt)
         }
         return true;
     }
-}
 
+    /// <summary>Verifies that every triangle has counter-clockwise winding.</summary>
+    public static bool VerifyWinding<T>(Triangulation<T> cdt)
+        where T : unmanaged, IFloatingPoint<T>, IMinMaxValue<T>, IRootFunctions<T>
+        => VerifyWinding(cdt.Vertices.Span, cdt.Triangles.Span);
+
+    /// <summary>Verifies that every triangle has counter-clockwise winding.</summary>
+    public static bool VerifyWinding<T>(
+        ReadOnlySpan<V2d<T>> vertices,
+        ReadOnlySpan<Triangle> triangles)
+        where T : unmanaged, IFloatingPoint<T>
+    {
+        foreach (ref readonly var triangle in triangles)
+        {
+            if ((uint)triangle.V0 >= (uint)vertices.Length ||
+                (uint)triangle.V1 >= (uint)vertices.Length ||
+                (uint)triangle.V2 >= (uint)vertices.Length ||
+                CdtUtils.Orient2D(
+                    vertices[triangle.V0],
+                    vertices[triangle.V1],
+                    vertices[triangle.V2]) <= T.Zero)
+                return false;
+        }
+        return true;
+    }
+}