Use unsafe and remove allocation in PossibleIntersection

Author: JimBobSquarePants

src/PolygonClipper/PolygonClipper.cs

@@ -5,6 +5,8 @@
 using System.Collections.Generic;
 using System.Diagnostics;
 using System.Diagnostics.CodeAnalysis;
+using System.Runtime.CompilerServices;
+using System.Runtime.InteropServices;
 
 namespace PolygonClipper;
 
@@ -162,6 +164,7 @@ public Polygon Run()
 
         SweepEvent? prevEvent;
         SweepEvent? nextEvent;
+        Span<SweepEvent> workspace = new SweepEvent[4];
         while (eventQueue.Count > 0)
         {
             SweepEvent sweepEvent = eventQueue.Dequeue();
@@ -188,7 +191,7 @@ public Polygon Run()
                 if (nextEvent != null)
                 {
                     // Check intersection with the next neighbor
-                    if (PossibleIntersection(sweepEvent, nextEvent, eventQueue) == 2)
+                    if (PossibleIntersection(sweepEvent, nextEvent, eventQueue, workspace) == 2)
                     {
                         ComputeFields(sweepEvent, prevEvent, operation);
                         ComputeFields(nextEvent, sweepEvent, operation);
@@ -199,7 +202,7 @@ public Polygon Run()
                 if (prevEvent != null)
                 {
                     // Check intersection with the previous neighbor
-                    if (PossibleIntersection(prevEvent, sweepEvent, eventQueue) == 2)
+                    if (PossibleIntersection(prevEvent, sweepEvent, eventQueue, workspace) == 2)
                     {
                         SweepEvent? prevPrevEvent = statusLine.Prev(prevEvent.PosSL);
                         ComputeFields(prevEvent, prevPrevEvent, operation);
@@ -218,7 +221,7 @@ public Polygon Run()
                 // Check intersection between neighbors
                 if (prevEvent != null && nextEvent != null)
                 {
-                    _ = PossibleIntersection(prevEvent, nextEvent, eventQueue);
+                    _ = PossibleIntersection(prevEvent, nextEvent, eventQueue, workspace);
                 }
 
                 statusLine.RemoveAt(it);
@@ -499,6 +502,10 @@ private static bool InResult(SweepEvent sweepEvent, BooleanOperation operation)
     /// <param name="le1">The first sweep event representing a line segment.</param>
     /// <param name="le2">The second sweep event representing a line segment.</param>
     /// <param name="eventQueue">The event queue to add new events to.</param>
+    /// <param name="workspace">
+    /// A scratch space for temporary storage of sweep events.
+    /// Must be at least 4 elements long to hold the events for the two segments and their associated other events.
+    /// </param>
     /// <returns>
     /// An integer indicating the result of the intersection:
     /// <list type="bullet">
@@ -514,7 +521,8 @@ private static bool InResult(SweepEvent sweepEvent, BooleanOperation operation)
     private static int PossibleIntersection(
         SweepEvent le1,
         SweepEvent le2,
-        StablePriorityQueue<SweepEvent, SweepEventComparer> eventQueue)
+        StablePriorityQueue<SweepEvent, SweepEventComparer> eventQueue,
+        Span<SweepEvent> workspace)
     {
         if (le1.OtherEvent == null || le2.OtherEvent == null)
         {
@@ -569,41 +577,55 @@ private static int PossibleIntersection(
         }
 
         // The line segments associated with le1 and le2 overlap.
-        // TODO: Rewrite this to avoid allocation.
-        List<SweepEvent> events = new(4);
         bool leftCoincide = le1.Point == le2.Point;
         bool rightCoincide = le1.OtherEvent.Point == le2.OtherEvent.Point;
 
-        // Populate the events
+        // Populate the events.
+        // The working buffer has a length of 4, which is sufficient to hold the events
+        // for the two segments and their associated other events.
+        // Events are assigned in a specific order to avoid overwriting shared references.
+        ref SweepEvent wRef = ref MemoryMarshal.GetReference(workspace);
         if (!leftCoincide)
         {
             if (comparer.Compare(le1, le2) > 0)
             {
-                events.Add(le2);
-                events.Add(le1);
+                Unsafe.Add(ref wRef, 0u) = le2;
+                Unsafe.Add(ref wRef, 1u) = le1;
             }
             else
             {
-                events.Add(le1);
-                events.Add(le2);
+                Unsafe.Add(ref wRef, 0u) = le1;
+                Unsafe.Add(ref wRef, 1u) = le2;
             }
-        }
 
-        if (!rightCoincide)
+            // Positions 0 and 1 contain the left events of the segments.
+            // Positions 2 and 3 will contain the right events of the segments.
+            if (!rightCoincide)
+            {
+                Unsafe.Add(ref wRef, 2u) = le1.OtherEvent;
+                Unsafe.Add(ref wRef, 3u) = le2.OtherEvent;
+            }
+            else
+            {
+                Unsafe.Add(ref wRef, 2u) = le2.OtherEvent;
+                Unsafe.Add(ref wRef, 3u) = le1.OtherEvent;
+            }
+        }
+        else if (leftCoincide && !rightCoincide)
         {
+            // Only the right endpoints differ, so we use positions 0 and 1 for their sorted order.
             if (comparer.Compare(le1.OtherEvent, le2.OtherEvent) > 0)
             {
-                events.Add(le2.OtherEvent);
-                events.Add(le1.OtherEvent);
+                Unsafe.Add(ref wRef, 0u) = le2.OtherEvent;
+                Unsafe.Add(ref wRef, 1u) = le1.OtherEvent;
             }
             else
             {
-                events.Add(le1.OtherEvent);
-                events.Add(le2.OtherEvent);
+                Unsafe.Add(ref wRef, 0u) = le1.OtherEvent;
+                Unsafe.Add(ref wRef, 1u) = le2.OtherEvent;
             }
         }
 
-        // Handle leftCoincide case
         if (leftCoincide)
         {
             le2.EdgeType = EdgeType.NonContributing;
@@ -613,30 +635,31 @@ private static int PossibleIntersection(
 
             if (leftCoincide && !rightCoincide)
             {
-                DivideSegment(events[1].OtherEvent, events[0].Point, eventQueue, comparer);
+                DivideSegment(Unsafe.Add(ref wRef, 1u).OtherEvent, Unsafe.Add(ref wRef, 0u).Point, eventQueue, comparer);
             }
 
             return 2;
         }
 
-        // Handle the rightCoincide case
         if (rightCoincide)
         {
-            DivideSegment(events[0], events[1].Point, eventQueue, comparer);
+            // Since leftCoincide is false, the first two workspace slots contain distinct left events.
+            DivideSegment(Unsafe.Add(ref wRef, 0u), Unsafe.Add(ref wRef, 1u).Point, eventQueue, comparer);
             return 3;
         }
 
         // Handle general overlapping case
-        if (events[0] != events[3].OtherEvent)
+        // At this point: workspace[0,1] = sorted left events, workspace[2,3] = sorted right events.
+        if (Unsafe.Add(ref wRef, 0u) != Unsafe.Add(ref wRef, 3u).OtherEvent)
         {
-            DivideSegment(events[0], events[1].Point, eventQueue, comparer);
-            DivideSegment(events[1], events[2].Point, eventQueue, comparer);
+            DivideSegment(Unsafe.Add(ref wRef, 0u), Unsafe.Add(ref wRef, 1u).Point, eventQueue, comparer);
+            DivideSegment(Unsafe.Add(ref wRef, 1u), Unsafe.Add(ref wRef, 2u).Point, eventQueue, comparer);
             return 3;
         }
 
         // One segment fully contains the other
-        DivideSegment(events[0], events[1].Point, eventQueue, comparer);
-        DivideSegment(events[3].OtherEvent, events[2].Point, eventQueue, comparer);
+        DivideSegment(Unsafe.Add(ref wRef, 0u), Unsafe.Add(ref wRef, 1u).Point, eventQueue, comparer);
+        DivideSegment(Unsafe.Add(ref wRef, 3u).OtherEvent, Unsafe.Add(ref wRef, 2u).Point, eventQueue, comparer);
         return 3;
     }
 
@@ -924,6 +947,7 @@ private static ReadOnlySpan<int> PrecomputeIterationOrder(List<SweepEvent> data)
         return map;
     }
 
+    [MethodImpl(MethodImplOptions.AggressiveInlining)]
     private static void MarkProcessed(SweepEvent sweepEvent, Span<bool> processed, int pos, int contourId)
     {
         processed[pos] = true;
@@ -994,11 +1018,9 @@ private static Contour InitializeContourFromContext(SweepEvent sweepEvent, Polyg
     /// starting from the given position.
     /// </summary>
     /// <param name="pos">The current position in the result events.</param>
-    /// <param name="resultEvents">The list of sweep events representing result segments.</param>
     /// <param name="processed">A list indicating whether each event at the corresponding index has been processed.</param>
-    /// <param name="originalPos">The original position to return if no unprocessed event is found.</param>
-    /// <param name="iterationMap"></param>
-    /// <param name="found"></param>
+    /// <param name="iterationMap">A precomputed map that indicates the next position to check for unprocessed events.</param>
+    /// <param name="found">A boolean indicating whether an unprocessed event was found.</param>
     /// <returns>The index of the next unprocessed position.</returns>
     /// <remarks>
     /// This method searches forward from the current position until it finds an unprocessed event with

tests/PolygonClipper.Benchmarks/ClippingLibraryComparison.cs

@@ -6,13 +6,23 @@
 using BenchmarkDotNet.Attributes;
 using Clipper2Lib;
 using GeoJSON.Text.Feature;
-using GeoJSON.Text.Geometry;
+using PolygonClipper.Tests;
 using PolygonClipper.Tests.TestCases;
 
-using GeoPolygon = GeoJSON.Text.Geometry.Polygon;
-
 namespace PolygonClipper.Benchmarks;
 
+/// <summary>
+/// <para>
+/// Benchmarks the performance of <c>PolygonClipper</c> against <c>Clipper2</c>
+/// for polygon union operations using a fixed input dataset.
+/// </para>
+/// <para>
+/// Note: Clipper2 produces an incorrect result for the test case used in this benchmark.
+/// As such, the Clipper2 timing is included for reference only and should not be considered
+/// a valid correctness comparison. This benchmark is intended to evaluate the performance
+/// of <c>PolygonClipper</c> in isolation.
+/// </para>
+/// </summary>
 public class ClippingLibraryComparison
 {
     private static readonly FeatureCollection Data = TestData.Generic.GetFeatureCollection("issue71.geojson");
@@ -25,17 +35,17 @@ public class ClippingLibraryComparison
     [GlobalSetup]
     public void Setup()
     {
-        (Polygon subject, Polygon clipping) = BuildPolygon();
+        (Polygon subject, Polygon clipping) = TestPolygonUtilities.BuildPolygon(Data);
         this.subject = subject;
         this.clipping = clipping;
 
-        (PathsD subject2, PathsD clipping2) = BuildPolygon2();
+        (PathsD subject2, PathsD clipping2) = TestPolygonUtilities.BuildClipper2Polygon(Data);
         this.subject2 = subject2;
         this.clipping2 = clipping2;
     }
 
     [Benchmark]
-    public Polygon Clipper() => PolygonClipper.Union(this.subject, this.clipping);
+    public Polygon PolygonClipper() => global::PolygonClipper.PolygonClipper.Union(this.subject, this.clipping);
 
     [Benchmark(Baseline = true)]
     public PathsD Clipper2()
@@ -47,108 +57,4 @@ public PathsD Clipper2()
         clipper2.Execute(ClipType.Union, FillRule.Positive, solution);
         return solution;
     }
-
-    public static (Polygon Subject, Polygon Clipping) BuildPolygon()
-    {
-        IGeometryObject subjectGeometry = Data.Features[0].Geometry;
-        IGeometryObject clippingGeometry = Data.Features[1].Geometry;
-
-        Polygon subject = ConvertToPolygon(subjectGeometry);
-        Polygon clipping = ConvertToPolygon(clippingGeometry);
-
-        return (subject, clipping);
-    }
-
-    public static (PathsD Subject, PathsD Clipping) BuildPolygon2()
-    {
-        IGeometryObject subjectGeometry = Data.Features[0].Geometry;
-        IGeometryObject clippingGeometry = Data.Features[1].Geometry;
-
-        PathsD subject = ConvertToPolygon2(subjectGeometry);
-        PathsD clipping = ConvertToPolygon2(clippingGeometry);
-
-        return (subject, clipping);
-    }
-
-    private static Polygon ConvertToPolygon(IGeometryObject geometry)
-    {
-        if (geometry is GeoPolygon geoJsonPolygon)
-        {
-            // Convert GeoJSON Polygon to our Polygon type
-            Polygon polygon = new();
-            foreach (LineString ring in geoJsonPolygon.Coordinates)
-            {
-                Contour contour = new();
-                foreach (IPosition xy in ring.Coordinates)
-                {
-                    contour.AddVertex(new Vertex(xy.Longitude, xy.Latitude));
-                }
-                polygon.Push(contour);
-            }
-
-            return polygon;
-        }
-        else if (geometry is MultiPolygon geoJsonMultiPolygon)
-        {
-            // Convert GeoJSON MultiPolygon to our Polygon type
-            Polygon polygon = new();
-            foreach (GeoPolygon geoPolygon in geoJsonMultiPolygon.Coordinates)
-            {
-                foreach (LineString ring in geoPolygon.Coordinates)
-                {
-                    Contour contour = new();
-                    foreach (IPosition xy in ring.Coordinates)
-                    {
-                        contour.AddVertex(new Vertex(xy.Longitude, xy.Latitude));
-                    }
-                    polygon.Push(contour);
-                }
-            }
-
-            return polygon;
-        }
-
-        throw new InvalidOperationException("Unsupported geometry type.");
-    }
-
-    private static PathsD ConvertToPolygon2(IGeometryObject geometry)
-    {
-        if (geometry is GeoPolygon geoJsonPolygon)
-        {
-            // Convert GeoJSON Polygon to our Polygon type
-            PathsD polygon = [];
-            foreach (LineString ring in geoJsonPolygon.Coordinates)
-            {
-                PathD contour = [];
-                foreach (IPosition xy in ring.Coordinates)
-                {
-                    contour.Add(new PointD(xy.Longitude, xy.Latitude));
-                }
-                polygon.Add(contour);
-            }
-
-            return polygon;
-        }
-        else if (geometry is MultiPolygon geoJsonMultiPolygon)
-        {
-            // Convert GeoJSON MultiPolygon to our Polygon type
-            PathsD polygon = [];
-            foreach (GeoPolygon geoPolygon in geoJsonMultiPolygon.Coordinates)
-            {
-                foreach (LineString ring in geoPolygon.Coordinates)
-                {
-                    PathD contour = [];
-                    foreach (IPosition xy in ring.Coordinates)
-                    {
-                        contour.Add(new PointD(xy.Longitude, xy.Latitude));
-                    }
-                    polygon.Add(contour);
-                }
-            }
-
-            return polygon;
-        }
-
-        throw new InvalidOperationException("Unsupported geometry type.");
-    }
 }

tests/PolygonClipper.Benchmarks/PolygonClipper.Benchmarks.csproj

@@ -31,11 +31,12 @@
   <ItemGroup>
     <Compile Include="..\PolygonClipper.Tests\TestData.cs" Link="TestData.cs" />
     <Compile Include="..\PolygonClipper.Tests\TestEnvironment.cs" Link="TestEnvironment.cs" />
+    <Compile Include="..\PolygonClipper.Tests\TestPolygonUtilities.cs" Link="TestPolygonUtilities.cs" />    
   </ItemGroup>
 
   <ItemGroup>
     <PackageReference Include="BenchmarkDotNet" Version="0.14.0" />
-    <PackageReference Include="Clipper2" Version="1.4.0" />
+    <PackageReference Include="Clipper2" Version="1.5.4" />
   </ItemGroup>
 
   <ItemGroup>

tests/PolygonClipper.Tests/GenericTestCases.cs

@@ -54,7 +54,7 @@ public void GenericTestCase(string testCaseFile)
         Polygon clipping = ConvertToPolygon(clippingGeometry);
 
 #pragma warning disable RCS1124 // Inline local variable
-        List<ExpectedResult> expectedResults = ExtractExpectedResults(data.Features.Skip(2).ToList(), data.Type);
+        List<ExpectedResult> expectedResults = ExtractExpectedResults([.. data.Features.Skip(2)], data.Type);
 #pragma warning restore RCS1124 // Inline local variable
 
         // ExpectedResult result = expectedResults[1];
@@ -157,7 +157,8 @@ private static List<ExpectedResult> ExtractExpectedResults(List<Feature> feature
 
             return new ExpectedResult
             {
-                Operation = operation, Coordinates = ConvertToPolygon(feature.Geometry as MultiPolygon)
+                Operation = operation,
+                Coordinates = ConvertToPolygon(feature.Geometry as MultiPolygon)
             };
         });
 

tests/PolygonClipper.Tests/PolygonClipper.Tests.csproj

@@ -23,6 +23,7 @@
     <PackageReference Include="xunit" Version="2.4.0" />
     <PackageReference Include="xunit.runner.visualstudio" Version="2.4.0" />
     <PackageReference Include="coverlet.collector" Version="1.2.0" />
+    <PackageReference Include="Clipper2" Version="1.5.4" />
   </ItemGroup>
 
   <ItemGroup>