diff --git a/LICENSE-predicates.txt b/LICENSE-predicates.txt new file mode 100644 index 0000000..fc4ab52 --- /dev/null +++ b/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. diff --git a/README.md b/README.md index 96166e2..912b4cd 100644 --- a/README.md +++ b/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 diff --git a/src/CDT.Core/CDT.Core.csproj b/src/CDT.Core/CDT.Core.csproj index 45001ec..f9d1eca 100644 --- a/src/CDT.Core/CDT.Core.csproj +++ b/src/CDT.Core/CDT.Core.csproj @@ -23,7 +23,7 @@ https://github.com/MichaCo/CDT.NET https://github.com/MichaCo/CDT.NET git - MPL-2.0 + MPL-2.0 AND BSD-3-Clause README.md true true @@ -33,7 +33,9 @@ + + diff --git a/src/CDT.Core/CdtUtils.cs b/src/CDT.Core/CdtUtils.cs index 8525ca8..1accb52 100644 --- a/src/CDT.Core/CdtUtils.cs +++ b/src/CDT.Core/CdtUtils.cs @@ -541,6 +541,70 @@ public static bool IsInCircumcircle(V2d p, V2d v1, V2d v2, V2d v3 Unsafe.As(ref Unsafe.AsRef(in p.X)), Unsafe.As(ref Unsafe.AsRef(in p.Y))) > 0f; } + /// + /// Returns whether lies strictly inside the circle + /// whose diameter is the segment from to . + /// + public static bool IsEncroachingOnEdge(V2d p, V2d a, V2d b) + where T : unmanaged, IFloatingPoint + => 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; + + /// Computes the area of triangle , , . + public static T Area(V2d a, V2d b, V2d c) + where T : unmanaged, IFloatingPoint + => T.Abs(Orient2D(a, b, c)) / (T.One + T.One); + + /// Computes the circumcenter of a non-degenerate triangle. + public static V2d Circumcenter(V2d a, V2d b, V2d c) + where T : unmanaged, IFloatingPoint + { + 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(a.X - c.X, a.Y - c.Y); + b = new V2d(b.X - c.X, b.Y - c.Y); + return new V2d( + c.X + (b.Y * aLenSq - a.Y * bLenSq) / denom, + c.Y + (a.X * bLenSq - b.X * aLenSq) / denom); + } + + /// Computes the smallest interior angle of a non-degenerate triangle, in radians. + public static T SmallestAngle(V2d a, V2d b, V2d c) + where T : unmanaged, IFloatingPoint, IRootFunctions + { + 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))); + } + + /// Converts an angle in degrees to radians. + public static T DegToRad(T degrees) + where T : IFloatingPoint + => degrees / T.CreateChecked(180) * T.Pi; + + internal static bool IsAngleAtApexSmall(V2d apex, V2d a, V2d b) + where T : unmanaged, IFloatingPoint, IRootFunctions + { + 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 a, T b, T t) + where T : IFloatingPoint + => (T.One - t) * a + t * b; + /// /// Computes the squared Euclidean distance between two points. /// diff --git a/src/CDT.Core/Predicates/PredicatesAdaptive.cs b/src/CDT.Core/Predicates/PredicatesAdaptive.cs index 153b709..a905253 100644 --- a/src/CDT.Core/Predicates/PredicatesAdaptive.cs +++ b/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); } + /// + /// Returns positive when point c lies strictly inside the circle with + /// diameter ab, zero when it lies on the circle, and negative when outside. + /// + 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 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 s1 = stackalloc double[4]; + int s1Len = TwoTwoDiff(acxtail, nbcx, acytail, bcy, s1); + Span s2 = stackalloc double[4]; + int s2Len = TwoTwoDiff(acx, nbcxtail, acy, bcytail, s2); + Span s3 = stackalloc double[4]; + int s3Len = TwoTwoDiff(acxtail, nbcxtail, acytail, bcytail, s3); + Span t1 = stackalloc double[8]; + int t1Len = ExpansionSum(b, bLen, s1, s1Len, t1); + Span t2 = stackalloc double[12]; + int t2Len = ExpansionSum(t1, t1Len, s2, s2Len, t2); + Span d = stackalloc double[16]; + int dLen = ExpansionSum(t2, t2Len, s3, s3Len, d); + return MostSignificant(d, dLen); + } + // ========================================================================= // Shewchuk/Lenthe floating-point expansion primitives (internal) // ========================================================================= diff --git a/src/CDT.Core/TopologyVerifier.cs b/src/CDT.Core/TopologyVerifier.cs index 355dcc7..ff986ea 100644 --- a/src/CDT.Core/TopologyVerifier.cs +++ b/src/CDT.Core/TopologyVerifier.cs @@ -53,5 +53,29 @@ public static bool VerifyTopology(Triangulation cdt) } return true; } -} + /// Verifies that every triangle has counter-clockwise winding. + public static bool VerifyWinding(Triangulation cdt) + where T : unmanaged, IFloatingPoint, IMinMaxValue, IRootFunctions + => VerifyWinding(cdt.Vertices.Span, cdt.Triangles.Span); + + /// Verifies that every triangle has counter-clockwise winding. + public static bool VerifyWinding( + ReadOnlySpan> vertices, + ReadOnlySpan triangles) + where T : unmanaged, IFloatingPoint + { + 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; + } +} diff --git a/src/CDT.Core/Triangulation.Refinement.cs b/src/CDT.Core/Triangulation.Refinement.cs new file mode 100644 index 0000000..64ed1f4 --- /dev/null +++ b/src/CDT.Core/Triangulation.Refinement.cs @@ -0,0 +1,454 @@ +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at https://mozilla.org/MPL/2.0/. + +using System.Numerics; + +namespace CDT; + +public sealed partial class Triangulation + where T : unmanaged, IFloatingPoint, IMinMaxValue, IRootFunctions +{ + /// + /// Refines bad-quality triangles by inserting Steiner vertices. + /// + /// Maximum number of vertices to add. + /// Criterion used to identify bad triangles. + /// Angle in radians or maximum area, depending on the criterion. + /// + /// Optional mutable set from a Collect* method. Its indices are updated + /// as refinement creates replacement triangles. + /// + /// Do not split edges or triangles at or below this length. + /// Counts of refinement attempts that could not be completed. + public Unrefined RefineTriangles( + int maxVerticesToInsert, + RefinementCriterion refinementCriterion = RefinementCriterion.SmallestAngle, + T? refinementThreshold = null, + HashSet? toEraseOrNull = null, + T? minEdgeLength = null) + { + if (IsFinalized) throw new TriangulationFinalizedException(); + ArgumentOutOfRangeException.ThrowIfNegative(maxVerticesToInsert); + + T threshold = refinementThreshold ?? CdtUtils.DegToRad(T.CreateChecked(20)); + T minLength = minEdgeLength ?? T.CreateChecked(1e-6); + if (_verticesCount > 0 && _kdTree == null) InitKdTree(); + + var unrefined = new Unrefined(); + int remainingVertexBudget = maxVerticesToInsert; + int steinerVerticesOffset = _verticesCount; + + ResolveEncroachedEdges( + new Queue(FindEncroachedFixedEdges()), + ref remainingVertexBudget, + steinerVerticesOffset, + circumcenter: null, + refinementCriterion, + threshold, + toEraseOrNull, + minLength, + ref unrefined); + + var badTriangles = new Queue(); + for (int iT = 0; iT < _trianglesCount; iT++) + { + var t = _triangles[iT]; + if (!CdtUtils.TouchesSuperTriangle(t) && + !(toEraseOrNull?.Contains(iT) ?? false) && + IsRefinementNeeded(t, refinementCriterion, threshold)) + badTriangles.Enqueue(iT); + } + + while (badTriangles.Count > 0 && remainingVertexBudget > 0) + { + int iT = badTriangles.Dequeue(); + if (toEraseOrNull?.Contains(iT) ?? false) continue; + + var t = _triangles[iT]; + if (!IsRefinementNeeded(t, refinementCriterion, threshold) || + (refinementCriterion == RefinementCriterion.SmallestAngle && IsSmallestAngleFixed(t))) + continue; + + int v0 = t.V0, v1 = t.V1, v2 = t.V2; + var p0 = _vertices[v0]; + var p1 = _vertices[v1]; + var p2 = _vertices[v2]; + T shortestEdge = T.Min( + CdtUtils.Distance(p0, p1), + T.Min(CdtUtils.Distance(p1, p2), CdtUtils.Distance(p2, p0))); + if (shortestEdge <= minLength) + { + unrefined.ShortEdgeTriangles++; + continue; + } + + if (CdtUtils.Orient2D(p0, p1, p2) == T.Zero) continue; + var circumcenter = CdtUtils.Circumcenter(p0, p1, p2); + int iCircumcenterTri = FindTriangleAt(circumcenter); + if (iCircumcenterTri == Indices.NoNeighbor) + { + unrefined.CircumcenterOutside++; + continue; + } + + var encroachedEdges = EdgesEncroachedBy(circumcenter, iCircumcenterTri); + int budgetBeforeSplits = remainingVertexBudget; + var newlyBadTriangles = ResolveEncroachedEdges( + new Queue(encroachedEdges), + ref remainingVertexBudget, + steinerVerticesOffset, + circumcenter, + refinementCriterion, + threshold, + toEraseOrNull, + minLength, + ref unrefined); + if (remainingVertexBudget == 0) break; + + if (remainingVertexBudget != budgetBeforeSplits || newlyBadTriangles.Count > 0) + { + foreach (int badTriangle in newlyBadTriangles) + badTriangles.Enqueue(badTriangle); + badTriangles.Enqueue(iT); + continue; + } + if (encroachedEdges.Count > 0) continue; + + var centerTriangle = _triangles[iCircumcenterTri]; + var location = CdtUtils.LocatePointTriangle( + circumcenter, + _vertices[centerTriangle.V0], + _vertices[centerTriangle.V1], + _vertices[centerTriangle.V2]); + if (location == PtTriLocation.OnVertex) + { + unrefined.CircumcenterOnVertex++; + continue; + } + + int iTopo = CdtUtils.IsOnEdge(location) + ? centerTriangle.GetNeighbor(CdtUtils.EdgeNeighborFromLocation(location)) + : Indices.NoNeighbor; + bool outsideInputHull = + CdtUtils.TouchesSuperTriangle(centerTriangle) && + (iTopo == Indices.NoNeighbor || + CdtUtils.TouchesSuperTriangle(_triangles[iTopo])); + if ((toEraseOrNull?.Contains(iCircumcenterTri) ?? false) || + (iTopo != Indices.NoNeighbor && (toEraseOrNull?.Contains(iTopo) ?? false)) || + outsideInputHull) + { + unrefined.CircumcenterOutside++; + continue; + } + + int iVert = _verticesCount; + AddNewVertex(circumcenter, Indices.NoNeighbor); + var stack = new Stack(4); + if (iTopo == Indices.NoNeighbor) + InsertVertexInsideTriangle(iVert, iCircumcenterTri, stack); + else + InsertVertexOnEdge(iVert, iCircumcenterTri, iTopo, handleFixedSplitEdge: true, stack); + EnsureDelaunayByEdgeFlips(iVert, stack); + TryAddVertexToLocator(iVert); + + int start = _vertTris[iVert]; + int curr = start; + do + { + var adjacent = _triangles[curr]; + if (!(toEraseOrNull?.Contains(curr) ?? false) && + !CdtUtils.TouchesSuperTriangle(adjacent) && + IsRefinementNeeded(adjacent, refinementCriterion, threshold)) + badTriangles.Enqueue(curr); + curr = adjacent.Next(iVert).triIndex; + } while (curr != start && curr != Indices.NoNeighbor); + remainingVertexBudget--; + } + + if (refinementCriterion == RefinementCriterion.SmallestAngle) + { + for (int iT = 0; iT < _trianglesCount; iT++) + { + var t = _triangles[iT]; + if (CdtUtils.TouchesSuperTriangle(t) || + (toEraseOrNull?.Contains(iT) ?? false)) + continue; + if (IsRefinementNeeded(t, refinementCriterion, threshold) && IsSmallestAngleFixed(t)) + unrefined.SharpFixedCorner++; + } + } + return unrefined; + } + + /// Finds fixed edges encroached by one of their opposed vertices. + public List FindEncroachedFixedEdges() + { + if (IsFinalized) throw new TriangulationFinalizedException(); + var encroached = new List(); + foreach (var edge in _fixedEdges) + { + if (IsEdgeEncroached(edge)) encroached.Add(edge); + } + encroached.Sort(static (a, b) => + { + int first = a.V1.CompareTo(b.V1); + return first != 0 ? first : a.V2.CompareTo(b.V2); + }); + return encroached; + } + + /// Finds triangles that do not meet the requested refinement criterion. + public List FindUnrefinedTriangles( + RefinementCriterion refinementCriterion = RefinementCriterion.SmallestAngle, + T? refinementThreshold = null) + { + T threshold = refinementThreshold ?? CdtUtils.DegToRad(T.CreateChecked(20)); + bool hasSuperTriangle = !IsFinalized; + var unrefined = new List(); + for (int iT = 0; iT < _trianglesCount; iT++) + { + var t = _triangles[iT]; + if (hasSuperTriangle && CdtUtils.TouchesSuperTriangle(t)) continue; + if (IsRefinementNeeded(t, refinementCriterion, threshold)) + unrefined.Add(iT); + } + return unrefined; + } + + private bool IsRefinementNeeded(Triangle t, RefinementCriterion criterion, T threshold) + { + var a = _vertices[t.V0]; + var b = _vertices[t.V1]; + var c = _vertices[t.V2]; + return criterion switch + { + RefinementCriterion.SmallestAngle => CdtUtils.SmallestAngle(a, b, c) < threshold, + RefinementCriterion.LargestArea => CdtUtils.Area(a, b, c) > threshold, + _ => throw new ArgumentOutOfRangeException(nameof(criterion)), + }; + } + + private bool IsSmallestAngleFixed(Triangle t) + { + int apexIndex = 0; + T shortestSide = CdtUtils.DistanceSquared(_vertices[t.V1], _vertices[t.V2]); + for (int i = 1; i < 3; i++) + { + T side = CdtUtils.DistanceSquared( + _vertices[t.GetVertex(CdtUtils.Ccw(i))], + _vertices[t.GetVertex(CdtUtils.Cw(i))]); + if (side < shortestSide) + { + shortestSide = side; + apexIndex = i; + } + } + int apex = t.GetVertex(apexIndex); + return _fixedEdges.Contains(new Edge(apex, t.GetVertex(CdtUtils.Ccw(apexIndex)))) && + _fixedEdges.Contains(new Edge(apex, t.GetVertex(CdtUtils.Cw(apexIndex)))); + } + + private bool IsEdgeEncroached(Edge edge) + { + var (iT, iTopo) = EdgeTriangles(edge.V1, edge.V2); + if (iT == Indices.NoNeighbor) return false; + var t = _triangles[iT]; + int edgeIndex = CdtUtils.EdgeNeighborIndex(t, edge.V1, edge.V2); + int opposed = t.GetVertex(CdtUtils.OpposedVertexIndex(edgeIndex)); + if (CdtUtils.IsEncroachingOnEdge(_vertices[opposed], _vertices[edge.V1], _vertices[edge.V2])) + return true; + return iTopo != Indices.NoNeighbor && + CdtUtils.IsEncroachingOnEdge( + _vertices[CdtUtils.OpposedVertex(_triangles[iTopo], iT)], + _vertices[edge.V1], + _vertices[edge.V2]); + } + + private List EdgesEncroachedBy(V2d vertex, int startTriangle) + { + var encroached = new HashSet(); + var traversed = new HashSet { startTriangle }; + var toTraverse = new Stack(); + toTraverse.Push(startTriangle); + while (toTraverse.Count > 0) + { + var t = _triangles[toTraverse.Pop()]; + for (int i = 0; i < 3; i++) + { + var edge = new Edge(t.GetVertex(CdtUtils.Ccw(i)), t.GetVertex(CdtUtils.Cw(i))); + if (_fixedEdges.Contains(edge)) + { + if (CdtUtils.IsEncroachingOnEdge(vertex, _vertices[edge.V1], _vertices[edge.V2])) + encroached.Add(edge); + continue; + } + + int iN = t.GetNeighbor(CdtUtils.OpposedNeighborIndex(i)); + if (iN == Indices.NoNeighbor || !traversed.Add(iN)) continue; + var neighbor = _triangles[iN]; + if (CdtUtils.IsInCircumcircle( + vertex, + _vertices[neighbor.V0], + _vertices[neighbor.V1], + _vertices[neighbor.V2])) + toTraverse.Push(iN); + } + } + var result = encroached.ToList(); + result.Sort(static (a, b) => + { + int first = a.V1.CompareTo(b.V1); + return first != 0 ? first : a.V2.CompareTo(b.V2); + }); + return result; + } + + private List ResolveEncroachedEdges( + Queue encroachedEdges, + ref int remainingVertexBudget, + int steinerVerticesOffset, + V2d? circumcenter, + RefinementCriterion refinementCriterion, + T badTriangleThreshold, + HashSet? toEraseOrNull, + T minEdgeLength, + ref Unrefined unrefined) + { + var badTriangles = new List(); + while (encroachedEdges.Count > 0 && remainingVertexBudget > 0) + { + var edge = encroachedEdges.Dequeue(); + if (!_fixedEdges.Contains(edge)) continue; + if (CdtUtils.Distance(_vertices[edge.V1], _vertices[edge.V2]) <= minEdgeLength) + { + unrefined.ShortEdges++; + continue; + } + + int? splitVertex = SplitEncroachedEdge( + edge, steinerVerticesOffset, toEraseOrNull, ref unrefined); + if (splitVertex is not int iSplitVert) continue; + remainingVertexBudget--; + + int start = _vertTris[iSplitVert]; + int iT = start; + do + { + var t = _triangles[iT]; + if (circumcenter.HasValue && + !(toEraseOrNull?.Contains(iT) ?? false) && + !CdtUtils.TouchesSuperTriangle(t) && + IsRefinementNeeded(t, refinementCriterion, badTriangleThreshold)) + badTriangles.Add(iT); + + for (int i = 0; i < 3; i++) + { + var triEdge = new Edge(t.GetVertex(i), t.GetVertex(CdtUtils.Cw(i))); + if (!_fixedEdges.Contains(triEdge)) continue; + if (IsEdgeEncroached(triEdge) || + (circumcenter.HasValue && + CdtUtils.IsEncroachingOnEdge( + circumcenter.Value, + _vertices[triEdge.V1], + _vertices[triEdge.V2]))) + encroachedEdges.Enqueue(triEdge); + } + iT = t.Next(iSplitVert).triIndex; + } while (iT != start && iT != Indices.NoNeighbor); + } + return badTriangles; + } + + private int? SplitEncroachedEdge( + Edge edge, + int steinerVerticesOffset, + HashSet? toEraseOrNull, + ref Unrefined unrefined) + { + var start = _vertices[edge.V1]; + var end = _vertices[edge.V2]; + var (iT, iTopo) = EdgeTriangles(edge.V1, edge.V2); + if (iT == Indices.NoNeighbor || iTopo == Indices.NoNeighbor) + { + unrefined.SplitVertexInvalid++; + return null; + } + + T split = T.One / (T.One + T.One); + if (edge.V1 < steinerVerticesOffset && + edge.V2 >= steinerVerticesOffset && + HasAnotherFixedEdgeAtSmallAngle(edge.V1, edge)) + { + T length = CdtUtils.Distance(start, end); + T halfLength = length / (T.One + T.One); + T nearestPowerOfTwo = T.One; + while (halfLength > nearestPowerOfTwo) + nearestPowerOfTwo += nearestPowerOfTwo; + while (halfLength < T.CreateChecked(0.75) * nearestPowerOfTwo) + nearestPowerOfTwo /= T.One + T.One; + split = nearestPowerOfTwo / length; + } + + var splitVertex = new V2d( + CdtUtils.Lerp(start.X, end.X, split), + CdtUtils.Lerp(start.Y, end.Y, split)); + if (!IsEdgeSplitVertexValid(splitVertex, iT, iTopo)) + { + unrefined.SplitVertexInvalid++; + return null; + } + + int iSplit = SplitFixedEdgeAt(edge, splitVertex, iT, iTopo); + if (toEraseOrNull != null) + { + if (toEraseOrNull.Contains(iT)) toEraseOrNull.Add(_trianglesCount - 2); + if (toEraseOrNull.Contains(iTopo)) toEraseOrNull.Add(_trianglesCount - 1); + } + return iSplit; + } + + private bool HasAnotherFixedEdgeAtSmallAngle(int vertex, Edge excludeEdge) + { + int other = excludeEdge.V1 == vertex ? excludeEdge.V2 : excludeEdge.V1; + int start = _vertTris[vertex]; + int iT = start; + do + { + var t = _triangles[iT]; + var (iNext, otherVertex) = t.Next(vertex); + var candidate = new Edge(vertex, otherVertex); + if (candidate != excludeEdge && + _fixedEdges.Contains(candidate) && + CdtUtils.IsAngleAtApexSmall( + _vertices[vertex], + _vertices[other], + _vertices[otherVertex])) + return true; + iT = iNext; + } while (iT != start && iT != Indices.NoNeighbor); + return false; + } + + private int FindTriangleAt(V2d position) + { + if (_trianglesCount == 0 || _verticesCount == 0) return Indices.NoNeighbor; + int startVertex = _kdTree?.Nearest(position.X, position.Y, _vertices) ?? 0; + int iT = WalkTriangles(startVertex, position); + if ((uint)iT < (uint)_trianglesCount) + { + var t = _triangles[iT]; + if (CdtUtils.LocatePointTriangle( + position, _vertices[t.V0], _vertices[t.V1], _vertices[t.V2]) != PtTriLocation.Outside) + return iT; + } + for (int i = 0; i < _trianglesCount; i++) + { + var t = _triangles[i]; + if (CdtUtils.LocatePointTriangle( + position, _vertices[t.V0], _vertices[t.V1], _vertices[t.V2]) != PtTriLocation.Outside) + return i; + } + return Indices.NoNeighbor; + } +} diff --git a/src/CDT.Core/Triangulation.cs b/src/CDT.Core/Triangulation.cs index 3437868..1fb7375 100644 --- a/src/CDT.Core/Triangulation.cs +++ b/src/CDT.Core/Triangulation.cs @@ -93,7 +93,7 @@ public InvalidEdgeSplitVertexException(Edge e1, Edge e2) /// Supports both constrained and conforming modes. /// /// Floating-point coordinate type (float or double). -public sealed class Triangulation +public sealed partial class Triangulation where T : unmanaged, IFloatingPoint, IMinMaxValue, IRootFunctions { // ------------------------------------------------------------------------- @@ -310,25 +310,14 @@ public void EraseSuperTriangle() { if (IsFinalized) throw new TriangulationFinalizedException(); if (_superGeomType != SuperGeometryType.SuperTriangle) return; - var toErase = new HashSet(); - for (int i = 0; i < _trianglesCount; i++) - { - if (CdtUtils.TouchesSuperTriangle(_triangles[i])) - toErase.Add(i); - } - FinalizeTriangulation(toErase); + FinalizeTriangulation(CollectSuperTriangle()); } /// Removes all outer triangles (flood-fill from super-triangle vertex until a constraint edge). public void EraseOuterTriangles() { if (IsFinalized) throw new TriangulationFinalizedException(); - if (_vertTris[0] == Indices.NoNeighbor) - throw new TriangulationException("No vertex triangle data – already finalized?"); - var seeds = new Stack(); - seeds.Push(_vertTris[0]); - var toErase = GrowToBoundary(seeds); - FinalizeTriangulation(toErase); + FinalizeTriangulation(CollectOuterTriangles()); } /// @@ -336,6 +325,37 @@ public void EraseOuterTriangles() /// using even-odd depth rule. /// public void EraseOuterTrianglesAndHoles() + { + if (IsFinalized) throw new TriangulationFinalizedException(); + FinalizeTriangulation(CollectOuterTrianglesAndHoles()); + } + + /// Collects triangles touching the super-triangle. + public HashSet CollectSuperTriangle() + { + if (IsFinalized) throw new TriangulationFinalizedException(); + var toErase = new HashSet(); + if (_superGeomType != SuperGeometryType.SuperTriangle) return toErase; + for (int i = 0; i < _trianglesCount; i++) + { + if (CdtUtils.TouchesSuperTriangle(_triangles[i])) toErase.Add(i); + } + return toErase; + } + + /// Collects triangles outside the constrained boundary. + public HashSet CollectOuterTriangles() + { + if (IsFinalized) throw new TriangulationFinalizedException(); + if (_vertTrisCount == 0 || _vertTris[0] == Indices.NoNeighbor) + throw new TriangulationException("No vertex triangle data – already finalized?"); + var seeds = new Stack(); + seeds.Push(_vertTris[0]); + return GrowToBoundary(seeds); + } + + /// Collects triangles outside the constrained boundary and in detected holes. + public HashSet CollectOuterTrianglesAndHoles() { if (IsFinalized) throw new TriangulationFinalizedException(); var depths = CalculateTriangleDepths(); @@ -344,7 +364,7 @@ public void EraseOuterTrianglesAndHoles() { if (depths[i] % 2 == 0) toErase.Add(i); } - FinalizeTriangulation(toErase); + return toErase; } /// @@ -593,7 +613,7 @@ private void InsertVertex_FlipFixedEdges(int iV, Stack stack, List fl out int itopo, out int iv2, out int iv3, out int iv4, out int n1, out int n2, out int n3, out int n4); - if (itopo != Indices.NoNeighbor && IsFlipNeeded(iV, iv2, iv3, iv4)) + if (itopo != Indices.NoNeighbor && IsFlipNeeded(iV, iv2, iv3, iv4, doFlipFixedEdges: true)) { var flippedEdge = new Edge(iv2, iv4); if (_fixedEdges.Contains(flippedEdge)) @@ -816,10 +836,14 @@ private void EdgeFlipInfo( } [MethodImpl(MethodImplOptions.AggressiveInlining)] - private bool IsFlipNeeded(int iV1, int iV2, int iV3, int iV4) + private bool IsFlipNeeded(int iV1, int iV2, int iV3, int iV4, bool doFlipFixedEdges = false) { - // Skip HashSet lookup when there are no fixed edges (pure vertex-insertion path). - if (_fixedEdges.Count > 0 && _fixedEdges.Contains(new Edge(iV2, iV4))) return false; + var oldEdge = new Edge(iV2, iV4); + if (!doFlipFixedEdges && _fixedEdges.Contains(oldEdge)) return false; + if (_fixedEdges.Count > 0 && + (IsSameOriginalEdge(new Edge(iV1, iV2), new Edge(iV2, iV3)) || + IsSameOriginalEdge(new Edge(iV3, iV4), new Edge(iV4, iV1)))) + return false; var v1 = _vertices[iV1]; var v2 = _vertices[iV2]; @@ -853,6 +877,18 @@ private bool IsFlipNeeded(int iV1, int iV2, int iV3, int iV4) return IsInCircumcircle(v1, v2, v3, v4); } + private bool IsSameOriginalEdge(Edge e1, Edge e2) + { + if (!_fixedEdges.Contains(e1) || !_fixedEdges.Contains(e2)) return false; + var original1 = _pieceToOriginals.TryGetValue(e1, out var pieces1) && pieces1.Count > 0 + ? pieces1[0] + : e1; + var original2 = _pieceToOriginals.TryGetValue(e2, out var pieces2) && pieces2.Count > 0 + ? pieces2[0] + : e2; + return original1 == original2; + } + private void FlipEdge( int iT, int iTopo, int v1, int v2, int v3, int v4, @@ -885,24 +921,30 @@ private void InsertEdgeIteration( int iA = edge.V1, iB = edge.V2; if (iA == iB) return; - if (HasEdge(iA, iB)) - { - FixEdge(edge, originalEdge); - return; - } - var a = _vertices[iA]; var b = _vertices[iB]; T distTol = _minDistToConstraintEdge == T.Zero ? T.Zero : _minDistToConstraintEdge * CdtUtils.Distance(a, b); + if (HasEdge(iA, iB)) + { + int iVSplit = OpposedVertexOnEdge(edge, distTol); + if (iVSplit != Indices.NoVertex) + { + remaining.Add(new Edge(iVSplit, iB)); + remaining.Add(new Edge(iA, iVSplit)); + return; + } + FixEdge(edge, originalEdge); + return; + } + var (iT, iVL, iVR) = IntersectedTriangle(iA, a, b, distTol); if (iT == Indices.NoNeighbor) { - var edgePart = new Edge(iA, iVL); - FixEdge(edgePart, originalEdge); remaining.Add(new Edge(iVL, iB)); + remaining.Add(new Edge(iA, iVL)); return; } @@ -1010,7 +1052,7 @@ private void HandleIntersectingEdgeStrategy( if (_fixedEdges.Contains(edgeLR)) { var newV = IntersectionPosition(_vertices[iA], _vertices[iB], _vertices[iVL], _vertices[iVR]); - if (!IsEdgeSplitVertexValid(newV, iT, iTopo, iVL, iVR)) + if (!IsEdgeSplitVertexValid(newV, iT, iTopo)) { throw new InvalidEdgeSplitVertexException( OriginalInputEdge(originalEdge), @@ -1038,8 +1080,21 @@ private void ConformToEdgeIteration( int iA = edge.V1, iB = edge.V2; if (iA == iB) return; + var a = _vertices[iA]; + var b = _vertices[iB]; + T distTol = _minDistToConstraintEdge == T.Zero + ? T.Zero + : _minDistToConstraintEdge * CdtUtils.Distance(a, b); + if (HasEdge(iA, iB)) { + int iVSplit = OpposedVertexOnEdge(edge, distTol); + if (iVSplit != Indices.NoVertex) + { + remaining.Add(new ConformToEdgeTask(new Edge(iVSplit, iB), originals, overlaps)); + remaining.Add(new ConformToEdgeTask(new Edge(iA, iVSplit), originals, overlaps)); + return; + } FixEdge(edge); if (overlaps > 0) _overlapCount[edge] = overlaps; if (originals.Count > 0 && edge != originals[0]) @@ -1047,20 +1102,11 @@ private void ConformToEdgeIteration( return; } - var a = _vertices[iA]; - var b = _vertices[iB]; - T distTol = _minDistToConstraintEdge == T.Zero - ? T.Zero - : _minDistToConstraintEdge * CdtUtils.Distance(a, b); - var (iT, iVleft, iVright) = IntersectedTriangle(iA, a, b, distTol); if (iT == Indices.NoNeighbor) { - var part = new Edge(iA, iVleft); - FixEdge(part); - if (overlaps > 0) _overlapCount[part] = overlaps; - InsertUnique(_pieceToOriginals.GetOrAdd(part), originals); remaining.Add(new ConformToEdgeTask(new Edge(iVleft, iB), originals, overlaps)); + remaining.Add(new ConformToEdgeTask(new Edge(iA, iVleft), originals, overlaps)); return; } @@ -1132,7 +1178,7 @@ private void HandleConformIntersecting( if (_fixedEdges.Contains(edgeLR)) { var newV = IntersectionPosition(_vertices[iA], _vertices[iB], _vertices[iVleft], _vertices[iVright]); - if (!IsEdgeSplitVertexValid(newV, iT, iTopo, iVleft, iVright)) + if (!IsEdgeSplitVertexValid(newV, iT, iTopo)) { throw new InvalidEdgeSplitVertexException( OriginalInputEdge(edge), @@ -1374,9 +1420,9 @@ private int SplitFixedEdgeAt(Edge edge, V2d splitVert, int iT, int iTopo) AddNewVertex(splitVert, Indices.NoNeighbor); var stack = new Stack(4); InsertVertexOnEdge(iSplit, iT, iTopo, handleFixedSplitEdge: false, stack); + SplitFixedEdge(edge, iSplit); TryAddVertexToLocator(iSplit); EnsureDelaunayByEdgeFlips(iSplit, stack); - SplitFixedEdge(edge, iSplit); return iSplit; } @@ -1393,31 +1439,20 @@ private int SplitFixedEdgeAt(Edge edge, V2d splitVert, int iT, int iTopo) /// selects the adjacent triangle it must be contained in; only the two other /// edges of that triangle are tested (robust predicates). /// - private bool IsEdgeSplitVertexValid(V2d splitVert, int iT, int iTopo, int iVL, int iVR) - { - // Orient the split edge as it appears (counter-clockwise) in iT. Locating - // the (floating-point-rounded) split vertex against it both tells whether - // the split is safe and which of the two triangles sharing the edge must - // contain the vertex: 'Left' is iT's interior side, 'Right' is iTopo's. - var tL = _triangles[iT]; - int sL = CdtUtils.EdgeNeighborIndex(tL, iVL, iVR); - var side = LocatePointLine(splitVert, _vertices[tL.GetVertex(sL)], _vertices[tL.GetVertex(CdtUtils.Ccw(sL))]); - if (side == PtLineLocation.OnLine) - { - return splitVert != _vertices[iVL] && splitVert != _vertices[iVR]; - } - var t = side == PtLineLocation.Left ? tL : _triangles[iTopo]; - - // The split vertex must not fall outside that triangle. Its relation to the - // split edge is already established by 'side', so only the two edges meeting - // at the apex (opposite the split edge) are tested. A point to the right of - // a counter-clockwise triangle's edge lies outside it. - int s = CdtUtils.EdgeNeighborIndex(t, iVL, iVR); - var from = _vertices[t.GetVertex(s)]; // split edge tail (CCW) - var to = _vertices[t.GetVertex(CdtUtils.Ccw(s))]; // split edge head (CCW) - var apex = _vertices[t.GetVertex(CdtUtils.Cw(s))]; // opposite the split edge - return LocatePointLine(splitVert, to, apex) != PtLineLocation.Right && - LocatePointLine(splitVert, apex, from) != PtLineLocation.Right; + private bool IsEdgeSplitVertexValid(V2d splitVert, int iT, int iTopo) + { + var t1 = _triangles[iT]; + int i = CdtUtils.OpposedVertexIndex(CdtUtils.NeighborIndex(t1, iTopo)); + var v1 = _vertices[t1.GetVertex(i)]; + var v2 = _vertices[t1.GetVertex(CdtUtils.Ccw(i))]; + var v4 = _vertices[t1.GetVertex(CdtUtils.Cw(i))]; + var t2 = _triangles[iTopo]; + var v3 = _vertices[t2.GetVertex( + CdtUtils.OpposedVertexIndex(CdtUtils.NeighborIndex(t2, iT)))]; + return LocatePointLine(splitVert, v1, v2) == PtLineLocation.Left && + LocatePointLine(splitVert, v2, v3) == PtLineLocation.Left && + LocatePointLine(splitVert, v3, v4) == PtLineLocation.Left && + LocatePointLine(splitVert, v4, v1) == PtLineLocation.Left; } /// @@ -1467,6 +1502,47 @@ private bool HasEdge(int va, int vb) return false; } + private (int iT, int iTopo) EdgeTriangles(int va, int vb) + { + int startTri = _vertTris[va]; + if (startTri == Indices.NoNeighbor) return (Indices.NoNeighbor, Indices.NoNeighbor); + + int iT = startTri; + do + { + var t = _triangles[iT]; + var (iNext, otherVertex) = t.Next(va); + if (otherVertex == vb) return (iT, iNext); + iT = iNext; + } while (iT != startTri && iT != Indices.NoNeighbor); + return (Indices.NoNeighbor, Indices.NoNeighbor); + } + + private int OpposedVertexOnEdge(Edge edge, T tolerance) + { + if (_fixedEdges.Contains(edge)) return Indices.NoVertex; + + var (iT, iTopo) = EdgeTriangles(edge.V1, edge.V2); + if (iT == Indices.NoNeighbor || iTopo == Indices.NoNeighbor) + return Indices.NoVertex; + + var a = _vertices[edge.V1]; + var b = _vertices[edge.V2]; + int[] opposedVertices = + [ + CdtUtils.OpposedVertex(_triangles[iT], iTopo), + CdtUtils.OpposedVertex(_triangles[iTopo], iT), + ]; + foreach (int iV in opposedVertices) + { + var v = _vertices[iV]; + if (LocatePointLine(v, a, b, tolerance) == PtLineLocation.OnLine && + CdtUtils.IsEncroachingOnEdge(v, a, b)) + return iV; + } + return Indices.NoVertex; + } + // ------------------------------------------------------------------------- // Finalization // ------------------------------------------------------------------------- @@ -1493,8 +1569,10 @@ private HashSet GrowToBoundary(Stack seeds) return traversed; } - private void FinalizeTriangulation(HashSet removedTriangles) + /// Removes the super-triangle and the specified triangles, then finalizes the mesh. + public void FinalizeTriangulation(IReadOnlySet removedTriangles) { + if (IsFinalized) throw new TriangulationFinalizedException(); _vertTrisCount = 0; if (_superGeomType == SuperGeometryType.SuperTriangle) @@ -1552,7 +1630,7 @@ private void RemapEdgesNoSuperTriangle(Dictionary> dict) foreach (var kv in updated) dict[kv.Key] = kv.Value; } - private void RemoveTriangles(HashSet removed) + private void RemoveTriangles(IReadOnlySet removed) { if (removed.Count == 0) return; // Build a flat bool[] for O(1) indexed lookup — replaces all HashSet.Contains calls diff --git a/src/CDT.Core/Types.cs b/src/CDT.Core/Types.cs index ba5119d..a2cfb0b 100644 --- a/src/CDT.Core/Types.cs +++ b/src/CDT.Core/Types.cs @@ -297,6 +297,38 @@ public enum IntersectingConstraintEdges DontCheck, } +/// Criterion used to identify triangles needing refinement. +public enum RefinementCriterion +{ + /// Refine triangles whose smallest angle is below the threshold. + SmallestAngle, + + /// Refine triangles whose area is above the threshold. + LargestArea, +} + +/// Counts of refinement attempts that could not be completed. +public struct Unrefined +{ + /// Triangles left unrefined because their shortest edge is too short. + public int ShortEdgeTriangles { get; internal set; } + + /// Triangles whose circumcenter lies outside the triangulated area. + public int CircumcenterOutside { get; internal set; } + + /// Triangles whose circumcenter coincides with an existing vertex. + public int CircumcenterOnVertex { get; internal set; } + + /// Triangles whose smallest angle is enclosed by fixed edges. + public int SharpFixedCorner { get; internal set; } + + /// Fixed edges left unsplit because their length is too short. + public int ShortEdges { get; internal set; } + + /// Fixed edges whose split vertex would break topology. + public int SplitVertexInvalid { get; internal set; } +} + /// Layer depth type (used for hole detection). public readonly struct LayerDepth { diff --git a/test/CDT.Tests/GroundTruthTests.cs b/test/CDT.Tests/GroundTruthTests.cs index 1f6876f..29d95f5 100644 --- a/test/CDT.Tests/GroundTruthTests.cs +++ b/test/CDT.Tests/GroundTruthTests.cs @@ -367,6 +367,7 @@ public void GroundTruth_Conforming(string inputFile, bool isTyped) cdt.EraseOuterTrianglesAndHoles(); Assert.True(TopologyVerifier.VerifyTopology(cdt)); + Assert.True(TopologyVerifier.VerifyWinding(cdt)); GroundTruthHelpers.AssertTopologyMatchesFile(cdt, expectedPath); } } diff --git a/test/CDT.Tests/PredicatesTests.cs b/test/CDT.Tests/PredicatesTests.cs index 84ea9b1..37146c4 100644 --- a/test/CDT.Tests/PredicatesTests.cs +++ b/test/CDT.Tests/PredicatesTests.cs @@ -46,6 +46,17 @@ public void InCircle_OutsideCircle_ReturnsNegative() double r = PredicatesAdaptive.InCircle(0, 0, 4, 0, 2, 4, 2, -2); Assert.True(r < 0); } + + [Theory] + [InlineData(0.5, 0.1, true)] + [InlineData(0.1, 0.5, false)] + [InlineData(0.5, 0.5, false)] + public void IsEncroachingOnEdge_UsesStrictDiametralCircle( + double x, double y, bool expected) + { + Assert.Equal(expected, CdtUtils.IsEncroachingOnEdge( + new V2d(x, y), new V2d(0, 0), new V2d(1, 0))); + } } /// Tests for triangle utils. @@ -144,4 +155,3 @@ public void EdgeSet_TreatsSymmetricEdgesAsSame() Assert.Contains(new Edge(2, 1), set); } } - diff --git a/test/CDT.Tests/RefinementTests.cs b/test/CDT.Tests/RefinementTests.cs new file mode 100644 index 0000000..0f4a725 --- /dev/null +++ b/test/CDT.Tests/RefinementTests.cs @@ -0,0 +1,171 @@ +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at https://mozilla.org/MPL/2.0/. + +namespace CDT.Tests; + +public sealed class RefinementTests +{ + [Fact] + public void RefineTriangles_SatisfiesSmallestAngleThreshold() + { + var cdt = CreateRectangle(100, 1); + var toErase = cdt.CollectOuterTrianglesAndHoles(); + TRefine(cdt, RefinementCriterion.SmallestAngle, CdtUtils.DegToRad(20.0), 10000, toErase); + + Assert.Empty(cdt.FindEncroachedFixedEdges()); + cdt.FinalizeTriangulation(toErase); + + Assert.True(TopologyVerifier.VerifyTopology(cdt)); + Assert.True(TopologyVerifier.VerifyWinding(cdt)); + Assert.True(cdt.Triangles.Length > 2); + Assert.Empty(cdt.FindUnrefinedTriangles( + RefinementCriterion.SmallestAngle, CdtUtils.DegToRad(20.0))); + } + + [Fact] + public void RefineTriangles_SatisfiesLargestAreaThreshold() + { + var cdt = CreateRectangle(10, 10); + var toErase = cdt.CollectOuterTrianglesAndHoles(); + TRefine(cdt, RefinementCriterion.LargestArea, 1.0, 1000, toErase); + + Assert.Empty(cdt.FindEncroachedFixedEdges()); + cdt.FinalizeTriangulation(toErase); + + Assert.True(TopologyVerifier.VerifyTopology(cdt)); + Assert.True(TopologyVerifier.VerifyWinding(cdt)); + Assert.Empty(cdt.FindUnrefinedTriangles(RefinementCriterion.LargestArea, 1.0)); + } + + [Fact] + public void RefineTriangles_ZeroAngleThresholdIsNoOp() + { + var cdt = CreateRectangle(10, 10); + int verticesBefore = cdt.Vertices.Length; + + var unrefined = cdt.RefineTriangles( + 100, RefinementCriterion.SmallestAngle, 0.0); + + Assert.Equal(verticesBefore, cdt.Vertices.Length); + Assert.Equal(default, unrefined); + } + + [Fact] + public void RefineTriangles_DoesNotExpandConvexHull() + { + var cdt = new Triangulation(); + cdt.InsertVertices([new(0, 0), new(10, 0), new(1, 1)]); + int verticesBefore = cdt.Vertices.Length; + + var unrefined = cdt.RefineTriangles( + 10, RefinementCriterion.SmallestAngle, CdtUtils.DegToRad(20.0)); + + Assert.Equal(verticesBefore, cdt.Vertices.Length); + Assert.True(unrefined.CircumcenterOutside > 0); + } + + [Fact] + public void RefineTriangles_RespectsVertexBudget() + { + var cdt = CreateRectangle(100, 1); + int verticesBefore = cdt.Vertices.Length; + + cdt.RefineTriangles( + 10, RefinementCriterion.SmallestAngle, CdtUtils.DegToRad(20.0)); + + Assert.InRange(cdt.Vertices.Length, verticesBefore, verticesBefore + 10); + Assert.True(TopologyVerifier.VerifyTopology(cdt)); + Assert.True(TopologyVerifier.VerifyWinding(cdt)); + } + + [Fact] + public void RefineTriangles_ReportsUnrefinableSharpFixedCorner() + { + double angle = CdtUtils.DegToRad(10.0); + var cdt = new Triangulation(); + cdt.InsertVertices( + [ + new V2d(0, 0), + new V2d(8, 0), + new V2d(4, 5), + new V2d(Math.Cos(angle), Math.Sin(angle)), + ]); + cdt.InsertEdges( + [ + new Edge(0, 1), new Edge(1, 2), new Edge(2, 3), new Edge(3, 0), + ]); + + var unrefined = cdt.RefineTriangles( + 100, RefinementCriterion.SmallestAngle, CdtUtils.DegToRad(20.0), + minEdgeLength: 0.05); + + Assert.Equal(1, unrefined.SharpFixedCorner); + Assert.True(TopologyVerifier.VerifyTopology(cdt)); + Assert.True(TopologyVerifier.VerifyWinding(cdt)); + } + + [Fact] + public void RefineTriangles_RejectsInvalidSplitVertexWithoutCorruptingTopology() + { + var cdt = new Triangulation( + VertexInsertionOrder.Auto, + IntersectingConstraintEdges.TryResolve, + 1e-8); + cdt.InsertVertices( + [ + new(7.6, 2.2), new(2, 1.2), new(3.4, 0.4), new(1.4, 1.4), + new(6.6, 6.2), new(3.8, 0.2), new(0.6, 5.8), new(2.8, 4), + new(3.6, 7.4), new(4.6, 4.8), new(1, 2.8), + ]); + cdt.InsertEdges( + [ + new(4, 7), new(0, 5), new(3, 6), new(1, 8), new(1, 5), + new(3, 5), new(2, 4), new(8, 10), new(0, 2), new(2, 5), + ]); + + var unrefined = cdt.RefineTriangles( + 300, RefinementCriterion.SmallestAngle, CdtUtils.DegToRad(25.0), + minEdgeLength: 1e-6); + + Assert.True(unrefined.SplitVertexInvalid > 0); + Assert.True(TopologyVerifier.VerifyTopology(cdt)); + Assert.True(TopologyVerifier.VerifyWinding(cdt)); + } + + [Fact] + public void CollectAndRefineMethods_RejectFinalizedTriangulation() + { + var cdt = CreateRectangle(10, 10); + cdt.EraseOuterTrianglesAndHoles(); + + Assert.Throws(() => cdt.CollectSuperTriangle()); + Assert.Throws(() => cdt.CollectOuterTriangles()); + Assert.Throws(() => cdt.CollectOuterTrianglesAndHoles()); + Assert.Throws(() => cdt.FinalizeTriangulation(new HashSet())); + Assert.Throws(() => cdt.FindEncroachedFixedEdges()); + Assert.Throws(() => cdt.RefineTriangles(100)); + } + + private static Triangulation CreateRectangle(double width, double height) + { + var cdt = new Triangulation(); + cdt.InsertVertices( + [ + new(0, 0), new(width, 0), new(width, height), new(0, height), + ]); + cdt.InsertEdges( + [ + new(0, 1), new(1, 2), new(2, 3), new(3, 0), + ]); + return cdt; + } + + private static void TRefine( + Triangulation cdt, + RefinementCriterion criterion, + double threshold, + int budget, + HashSet toErase) + => cdt.RefineTriangles(budget, criterion, threshold, toErase); +} diff --git a/test/CDT.Tests/TriangulationTests.cs b/test/CDT.Tests/TriangulationTests.cs index d3a6a41..9516a42 100644 --- a/test/CDT.Tests/TriangulationTests.cs +++ b/test/CDT.Tests/TriangulationTests.cs @@ -30,6 +30,7 @@ public void InsertFourPoints_TopologicallyValid() cdt.InsertVertices([Pt(0, 0), Pt(1, 1), Pt(3, 1), Pt(3, 0)]); Assert.True(TopologyVerifier.VerifyTopology(cdt)); + Assert.True(TopologyVerifier.VerifyWinding(cdt)); // 4 user vertices + 3 super-triangle vertices Assert.Equal(7, cdt.Vertices.Length); Assert.Empty(cdt.FixedEdges); @@ -37,6 +38,15 @@ public void InsertFourPoints_TopologicallyValid() Assert.Equal(9, cdt.Triangles.Length); } + [Fact] + public void VerifyWinding_RejectsClockwiseTriangle() + { + var vertices = new[] { Pt(0, 0), Pt(1, 0), Pt(0, 1) }; + var triangle = new[] { new Triangle(0, 2, 1, -1, -1, -1) }; + + Assert.False(TopologyVerifier.VerifyWinding(vertices, triangle)); + } + [Fact] public void EraseSuperTriangle_LeavesConvexHull() { @@ -263,6 +273,59 @@ public void ConformToEdges_AfterFinalized_ThrowsTriangulationFinalizedException( cdt.ConformToEdges([new Edge(0, 2)])); } + [Theory] + [InlineData(false)] + [InlineData(true)] + public void ConstraintEdge_NearExistingVertex_SplitsAtVertex(bool conforming) + { + var cdt = new Triangulation( + VertexInsertionOrder.AsProvided, + IntersectingConstraintEdges.NotAllowed, + F(1e-3)); + cdt.InsertVertices( + [ + Pt(1, 0), Pt(0, 1), Pt(0, 0), + Pt(0.23263432516898214, 0.76736567483101781), + Pt(0.13422178907162741, 0), + ]); + var edges = new[] + { + new Edge(0, 1), new Edge(1, 2), new Edge(2, 0), + new Edge(3, 4), new Edge(3, 4), + }; + + if (conforming) cdt.ConformToEdges(edges); + else cdt.InsertEdges(edges); + cdt.EraseOuterTrianglesAndHoles(); + + Assert.Equal(3, cdt.Triangles.Length); + Assert.Contains(new Edge(0, 3), cdt.FixedEdges); + Assert.Contains(new Edge(1, 3), cdt.FixedEdges); + } + + [Theory] + [InlineData(false)] + [InlineData(true)] + public void ConstraintEdge_NearVertexOutsideDiametralCircle_IsNotSplit(bool conforming) + { + var cdt = new Triangulation( + VertexInsertionOrder.AsProvided, + IntersectingConstraintEdges.NotAllowed, + F(1e-3)); + cdt.InsertVertices( + [ + Pt(0, 0), Pt(1, 0), Pt(1e-7, 5e-4), Pt(0.5, 1), Pt(-1, 1), + ]); + var edges = new[] { new Edge(0, 1) }; + + if (conforming) cdt.ConformToEdges(edges); + else cdt.InsertEdges(edges); + cdt.EraseSuperTriangle(); + + Assert.Single(cdt.FixedEdges); + Assert.Contains(new Edge(0, 1), cdt.FixedEdges); + } + [Fact] public void EraseSuperTriangle_AfterFinalized_ThrowsTriangulationFinalizedException() { diff --git a/test/CDT.Tests/expected/issue-142-double-hanging-edge__conforming_auto_ignore_auto.txt b/test/CDT.Tests/expected/issue-142-double-hanging-edge__conforming_auto_ignore_auto.txt index a79af39..1d4af85 100644 --- a/test/CDT.Tests/expected/issue-142-double-hanging-edge__conforming_auto_ignore_auto.txt +++ b/test/CDT.Tests/expected/issue-142-double-hanging-edge__conforming_auto_ignore_auto.txt @@ -1,64 +1,89 @@ -18 -0 7 9 5 16 1 +26 +0 7 9 5 20 1 0 9 12 0 4294967295 4 -0 10 1 3 4294967295 4294967295 +0 10 24 3 23 4294967295 0 11 10 4 4294967295 2 0 12 11 1 4294967295 3 -0 19 7 4294967295 17 0 -2 5 15 7 13 4294967295 -2 6 5 8 4294967295 6 -2 7 6 9 4294967295 7 -2 13 7 4294967295 15 8 -3 18 4 4294967295 12 4294967295 -4 17 5 12 14 4294967295 -4 18 17 10 4294967295 11 -5 16 15 14 4294967295 6 -5 17 16 11 4294967295 13 -7 13 14 9 4294967295 16 -7 14 9 15 4294967295 0 -7 19 8 5 4294967295 4294967295 +0 23 7 4294967295 22 0 +1 27 25 4294967295 24 4294967295 +2 5 15 8 16 4294967295 +2 6 5 9 4294967295 7 +2 7 6 10 4294967295 8 +2 13 7 4294967295 17 9 +3 20 4 4294967295 13 4294967295 +4 18 5 13 15 4294967295 +4 20 18 11 4294967295 12 +5 17 21 15 4294967295 16 +5 18 17 12 4294967295 14 +5 21 15 14 4294967295 7 +7 13 16 10 4294967295 19 +7 14 19 21 4294967295 20 +7 16 22 17 4294967295 21 +7 19 9 18 4294967295 0 +7 22 14 19 4294967295 18 +7 23 8 5 4294967295 4294967295 +10 26 24 4294967295 25 2 +24 25 27 25 6 4294967295 +24 26 25 23 4294967295 24 -20 -0 1 -0 19 -1 10 +28 +0 23 +0 24 +1 25 +1 27 2 13 2 15 3 4 -3 18 +3 20 4 5 5 6 6 7 7 8 -8 19 +8 23 9 12 -9 14 +9 19 10 11 +10 26 11 12 -13 14 -15 16 -16 17 +13 16 +14 19 +14 22 +15 21 +16 22 17 18 +17 21 +18 20 +24 27 +25 26 0 -14 -0 19 +27 +0 23 1 0 8 +0 24 + 1 + 0 1 1 10 1 1 2 +1 25 + 1 + 1 2 +1 27 + 1 + 0 1 2 13 1 1 2 2 15 1 2 3 -3 18 +3 20 1 2 3 -8 19 +8 23 1 0 8 9 12 @@ -67,21 +92,51 @@ 9 14 1 1 2 +9 19 + 1 + 1 2 10 11 1 1 2 +10 26 + 1 + 1 2 11 12 1 1 2 13 14 1 1 2 -15 16 +13 16 1 - 2 3 -16 17 + 1 2 +14 16 + 1 + 1 2 +14 19 + 1 + 1 2 +14 22 + 1 + 1 2 +15 21 1 2 3 +16 22 + 1 + 1 2 17 18 1 2 3 +17 21 + 1 + 2 3 +18 20 + 1 + 2 3 +24 27 + 1 + 0 1 +25 26 + 1 + 1 2