Use lazy iterator to ensure all rings in create_[a/b]_list are closed

src/GeometryOps.jl

@@ -21,6 +21,7 @@ const Edge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T
 
 include("types.jl")
 include("primitives.jl")
+include("lazy_wrappers.jl")
 include("utils.jl")
 include("not_implemented_yet.jl")
 

src/lazy_wrappers.jl

@@ -0,0 +1,119 @@
+#=
+# Lazy wrappers
+
+These wrappers lazily apply some fixes like closing rings.
+=#
+
+abstract type AbstractLazyWrapper{GeomType} end
+
+struct LazyClosedRing{GeomType} <: AbstractLazyWrapper{GeomType}
+    ring::GeomType
+    function LazyClosedRing(ring)
+        LazyClosedRing(GI.trait(ring), ring)
+    end
+    function LazyClosedRing{GeomType}(ring::GeomType) where GeomType
+        new{GeomType}(ring)
+    end
+end
+
+
+
+function LazyClosedRing(::GI.AbstractCurveTrait, ring::GeomType) where GeomType
+    LazyClosedRing{GeomType}(ring)
+end
+
+# GeoInterface implementation
+GI.geomtrait(::LazyClosedRing) = GI.LinearRingTrait()
+GI.is3d(wrapper::LazyClosedRing) = GI.is3d(wrapper.ring)
+GI.ismeasured(wrapper::LazyClosedRing) = GI.ismeasured(wrapper.ring)
+GI.isclosed(::LazyClosedRing) = true
+
+function GI.npoint(wrapper::LazyClosedRing)
+    ring_npoints = GI.npoint(wrapper.ring)
+    if GI.getpoint(wrapper.ring, 1) == GI.getpoint(wrapper.ring, ring_npoints)
+        return ring_npoints
+    else
+        return ring_npoints + 1 # account for closing
+    end
+end
+
+function GI.getpoint(wrapper::LazyClosedRing)
+    return LazyClosedRingTuplePointIterator(wrapper)
+end
+
+function GI.getpoint(wrapper::LazyClosedRing, i::Integer)
+    ring_npoint = GI.npoint(wrapper.ring)
+    if i ≤ ring_npoint
+        return GI.getpoint(wrapper.ring, i)
+    elseif i == ring_npoint + 1
+        if GI.getpoint(wrapper.ring, 1) == GI.getpoint(wrapper.ring, ring_npoint)
+            return throw(BoundsError(wrapper.ring, i))
+        else
+            return GI.getpoint(wrapper.ring, 1)
+        end
+    else
+        return throw(BoundsError(wrapper.ring, i))
+    end
+end
+
+function tuples(wrapper::LazyClosedRing)
+    return collect(LazyClosedRingTuplePointIterator(wrapper))
+end
+
+struct LazyClosedRingTuplePointIterator{hasZ, hasM, GeomType}
+    ring::GeomType
+    closed::Bool
+end
+
+function LazyClosedRingTuplePointIterator(ring::LazyClosedRing)
+    geom = ring.ring
+    if GI.isempty(geom)
+        return LazyClosedRingTuplePointIterator{GI.is3d(geom), GI.ismeasured(geom), typeof(geom)}(geom, true)
+    end
+    isclosed = GI.getpoint(geom, 1) == GI.getpoint(geom, GI.npoint(geom))
+    return LazyClosedRingTuplePointIterator{GI.is3d(geom), GI.ismeasured(geom), typeof(geom)}(geom, isclosed)
+end
+
+# Base iterator interface
+Base.IteratorSize(::LazyClosedRingTuplePointIterator) = Base.HasLength()
+Base.length(iter::LazyClosedRingTuplePointIterator) = GI.npoint(iter.ring) + !iter.closed
+Base.IteratorEltype(::LazyClosedRingTuplePointIterator) = Base.HasEltype()
+function Base.eltype(::LazyClosedRingTuplePointIterator{hasZ, hasM}) where {hasZ, hasM}
+    if !hasZ && !hasM
+        Tuple{Float64, Float64}
+    elseif hasZ ⊻ hasM
+        Tuple{Float64, Float64, Float64}
+    else # hasZ && hasM
+        Tuple{Float64, Float64, Float64, Float64}
+    end
+end
+
+function Base.iterate(iter::LazyClosedRingTuplePointIterator)
+    if GI.isempty(iter.ring)
+        return nothing
+    else
+        return (GI.getpoint(iter.ring, 1), 1)
+    end
+end
+
+function Base.iterate(iter::LazyClosedRingTuplePointIterator, state)
+    ring_npoint = GI.npoint(iter.ring)
+    if iter.closed 
+        if state == ring_npoint
+            return nothing
+        else
+            return (GI.getpoint(iter.ring, state + 1), state + 1)
+        end
+    else
+        if state < ring_npoint
+            return (GI.getpoint(iter.ring, state + 1), state + 1)
+        elseif state == ring_npoint 
+            return (GI.getpoint(iter.ring, 1), state + 1)
+        elseif state == ring_npoint + 1
+            return nothing
+        else
+            throw(BoundsError(iter.ring, state))
+        end
+    end
+end
+

src/methods/clipping/difference.jl

@@ -51,8 +51,8 @@ function _difference(
     exact, kwargs...
 ) where T
     # Get the exterior of the polygons
-    ext_a = GI.getexterior(poly_a)
-    ext_b = GI.getexterior(poly_b)
+    ext_a = LazyClosedRing(GI.getexterior(poly_a))
+    ext_b = LazyClosedRing(GI.getexterior(poly_b))
     # Find the difference of the exterior of the polygons
     a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b, _diff_delay_cross_f, _diff_delay_bounce_f; exact)
     polys = _trace_polynodes(T, a_list, b_list, a_idx_list, _diff_step, poly_a, poly_b)

src/methods/clipping/intersection.jl

@@ -67,8 +67,8 @@ function _intersection(
     exact, kwargs...,
 ) where {T}
     # First we get the exteriors of 'poly_a' and 'poly_b'
-    ext_a = GI.getexterior(poly_a)
-    ext_b = GI.getexterior(poly_b)
+    ext_a = LazyClosedRing(GI.getexterior(poly_a))
+    ext_b = LazyClosedRing(GI.getexterior(poly_b))
     # Then we find the intersection of the exteriors
     a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b, _inter_delay_cross_f, _inter_delay_bounce_f; exact)
     polys = _trace_polynodes(T, a_list, b_list, a_idx_list, _inter_step, poly_a, poly_b)

src/methods/clipping/union.jl

@@ -51,8 +51,8 @@ function _union(
     exact, kwargs...,
 ) where T
     # First, I get the exteriors of the two polygons
-    ext_a = GI.getexterior(poly_a)
-    ext_b = GI.getexterior(poly_b)
+    ext_a = LazyClosedRing(GI.getexterior(poly_a))
+    ext_b = LazyClosedRing(GI.getexterior(poly_b))
     # Then, I get the union of the exteriors
     a_list, b_list, a_idx_list = _build_ab_list(T, ext_a, ext_b, _union_delay_cross_f, _union_delay_bounce_f; exact)
     polys = _trace_polynodes(T, a_list, b_list, a_idx_list, _union_step, poly_a, poly_b)

test/lazy_wrappers.jl

@@ -0,0 +1,95 @@
+# Tests for lazy wrappers
+
+# - Test that return type is inferred
+# - Test that results are correct
+# - Test that lazy evaluation is performed
+# - Test for proper error handling
+# - Test compatibility with different input types
+
+using Test
+using GeometryOps
+import GeoInterface as GI, GeometryOps as GO
+using ..TestHelpers
+
+@testset "LazyClosedRing" begin
+    # Helper function to create a simple LineString
+    create_linestring(closed=false) = closed ? GI.LineString([
+        (0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0),
+
+    ]) : GI.LineString([
+        (0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)
+    ])
+
+    @testset "Type inference" begin
+        ls = create_linestring()
+        wrapped = GO.LazyClosedRing(ls)
+        @inferred GI.npoint(wrapped)
+        @inferred GI.getpoint(wrapped, 1)
+        @inferred collect(GI.getpoint(wrapped))
+    end
+
+    @testset "Correctness" begin
+        ls = create_linestring()
+        wrapped = GO.LazyClosedRing(ls)
+
+        @test GI.npoint(wrapped) == 5
+        @test GI.getpoint(wrapped, 1) == (0.0, 0.0)
+        @test GI.getpoint(wrapped, 5) == (0.0, 0.0)
+        @test collect(GI.getpoint(wrapped)) == [
+            (0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)
+        ]
+    end
+
+    @testset "Lazy evaluation" begin
+        ls = create_linestring()
+        wrapped = GO.LazyClosedRing(ls)
+
+        # Check that the original linestring is not modified
+        @test GI.npoint(ls) == 4
+        @test GI.npoint(wrapped) == 5
+    end
+
+    @testset "Error handling" begin
+        ls = create_linestring()
+        wrapped = GO.LazyClosedRing(ls)
+
+        @test_throws BoundsError GI.getpoint(wrapped, 0)
+        @test_throws BoundsError GI.getpoint(wrapped, 6)
+    end
+
+    @testset "Compatibility with different input types" begin
+        # Test with a closed LineString
+        closed_ls = create_linestring(true)
+        wrapped_closed = GO.LazyClosedRing(closed_ls)
+        @test GI.npoint(wrapped_closed) == 5
+        @test collect(GI.getpoint(wrapped_closed)) == GI.getpoint(closed_ls)
+
+        # Test with a 3D LineString
+        ls_3d = GI.LineString([(x, y, 0.0) for (x, y) in GI.getpoint(create_linestring())])
+        wrapped_3d = GO.LazyClosedRing(ls_3d)
+        @test GI.is3d(wrapped_3d) == true
+        @test GI.npoint(wrapped_3d) == 5
+        @test GI.getpoint(wrapped_3d, 5) == (0.0, 0.0, 0.0)
+
+        # Test with a measured LineString
+        ls_measured = GI.LineString([GI.Point{false, true}(x, y, 0.0) for (x, y) in GI.getpoint(create_linestring())])
+        wrapped_measured = GO.LazyClosedRing(ls_measured)
+        @test GI.ismeasured(wrapped_measured) == true
+        @test GI.npoint(wrapped_measured) == 5
+    end
+end
+
+@testset "LazyClosedRingTuplePointIterator" begin
+    # Helper function to create a simple LineString
+    create_linestring(closed=false) = closed ? GI.LineString([
+        (0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0),
+
+    ]) : GI.LineString([
+        (0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)
+    ])
+    @testset "Type inference" begin
+        ls = create_linestring()
+        wrapped = GO.LazyClosedRing(ls)
+        @inferred Tuple{Float64, Float64} eltype(GO.LazyClosedRingTuplePointIterator(wrapped))
+    end
+end

test/methods/clipping/polygon_clipping.jl

@@ -106,6 +106,8 @@ p54 = GI.Polygon([[(2.5, 2.5), (2.5, 7.5), (7.5, 7.5), (7.5, 2.5), (2.5, 2.5)],
 p55 = GI.Polygon([[(5.0, 0.25), (5.0, 5.0), (9.5, 5.0), (9.5, 0.25), (5.0, 0.25)],
     [(6.0, 3.0), (6.0, 4.0), (7.0, 4.0), (7.0, 3.0), (6.0, 3.0)],
     [(7.5, 0.5), (7.5, 2.5), (9.25, 2.5), (9.25, 0.5), (7.5, 0.5)]])
+p56 = GI.Polygon([[(0.0, 0.0), (0.0, 10.0), (10.0, 10.0), (10.0, 0.0)]]) # polygons with unclosed rings
+p57 = GI.Polygon([[(0.0, 0.0), (0.0, 11.0), (11.0, 11.0), (11.0, 0.0)]])
 
 mp1 = GI.MultiPolygon([p1])
 mp2 = GI.MultiPolygon([p2])
@@ -160,7 +162,8 @@ test_pairs = [
     (p34, mp3, "p34", "mp3", "Polygon overlaps with multipolygon, where on of the sub-polygons is equivalent"),
     (mp3, p35, "mp3", "p35", "Mulitpolygon where just one sub-polygon touches other polygon"),
     (p35, mp3, "p35", "mp3", "Polygon that touches just one of the sub-polygons of multipolygon"),
-    (mp4, mp3, "mp4", "mp3", "Two multipolygons, which with two sub-polygons, where some of the sub-polygons intersect and some don't")
+    (mp4, mp3, "mp4", "mp3", "Two multipolygons, which with two sub-polygons, where some of the sub-polygons intersect and some don't"),
+    # (p56, p57, "p56", "p57", "Polygons with unclosed rings (#191)"), # TODO: `difference` doesn't work yet on p56 and p57 in that order (#193)
 ]
 
 const ϵ = 1e-10
@@ -215,3 +218,33 @@ end
 @testset "Intersection" begin test_clipping(GO.intersection, LG.intersection, "intersection") end
 @testset "Union" begin test_clipping(GO.union, LG.union, "union") end
 @testset "Difference" begin test_clipping(GO.difference, LG.difference, "difference") end
+
+
+@testset "Lazy closed ring enumerator" begin
+    # first test an open ring
+    p = GI.Polygon([[(0.0, 0.0), (0.0, 10.0), (10.0, 10.0), (10.0, 0.0)]])
+    cl = collect(GO.LazyClosedRingTuplePointIterator(GO.LazyClosedRing(GI.getring(p, 1))))
+    @test length(cl) == 5
+    @test cl[end][2] == cl[1][2]
+    # then test a closed ring
+    p2 = GI.Polygon([[(0.0, 0.0), (0.0, 10.0), (10.0, 10.0), (10.0, 0.0), (0.0, 0.0)]])
+    cl2 = collect(GO.LazyClosedRingTuplePointIterator(GO.LazyClosedRing(GI.getring(p2, 1))))
+    @test length(cl2) == 5
+    @test cl2[end][2] == cl2[1][2]
+    @test all(cl .== cl2)
+    # TODO: `difference` doesn't work yet on p56 and p57 in that order,
+    # so we do some tests here
+    p_intersection = only(GO.intersection(p56, p57; target = GI.PolygonTrait()))
+    p_union = only(GO.union(p56, p57; target = GI.PolygonTrait()))
+    @test GI.extent(p_intersection) == GI.extent(p56)
+    @test GI.extent(p_union) == GI.extent(p57)
+    @test GO.equals(p_intersection, p56)
+    @test GO.equals(p_union, p57)
+    # Notice how the order of polygons is different here
+    p_diff = only(GO.difference(p57, p56; target = GI.PolygonTrait()))
+    @test GI.extent(p_diff) == GI.extent(p57)
+    p_lg_diff = LG.difference(GO.fix(p57), GO.fix(p56))
+    # The point orders differ so we have to run GO intersection again.
+    # We could test area of difference also like `compare_GO_LG_clipping` does.
+    @test GO.equals(p_diff, only(GO.intersection(p_diff, p_lg_diff; target = GI.PolygonTrait())))
+end

test/runtests.jl

@@ -5,6 +5,7 @@ using SafeTestsets
 include("helpers.jl")
 
 @safetestset "Primitives" begin include("primitives.jl") end
+@safetestset "Lazy Wrappers" begin include("lazy_wrappers.jl") end
 # Methods
 @safetestset "Angles" begin include("methods/angles.jl") end
 @safetestset "Area" begin include("methods/area.jl") end