From e39e48facb571adacf0ea214258ecd0056f5ef2e Mon Sep 17 00:00:00 2001
From: Nemanja Jovanovic <nemanja.jovanovic@proton.me>
Date: Mon, 12 Aug 2024 18:57:33 +0300
Subject: [PATCH 1/6] Add Weiler-Atherton clipping algorithm

---
 docs/src/algorithms/clipping.md |  24 ++
 src/Meshes.jl                   |   1 +
 src/clipping.jl                 |   1 +
 src/clipping/weileratherton.jl  | 432 ++++++++++++++++++++++++++++++++
 test/clipping.jl                | 173 +++++++++++++
 5 files changed, 631 insertions(+)
 create mode 100644 src/clipping/weileratherton.jl

diff --git a/docs/src/algorithms/clipping.md b/docs/src/algorithms/clipping.md
index c2655598f..947126b58 100644
--- a/docs/src/algorithms/clipping.md
+++ b/docs/src/algorithms/clipping.md
@@ -28,6 +28,30 @@ other = Box((0,1), (3,7))
 # clipped polygon
 clipped = clip(poly, other, SutherlandHodgmanClipping())
 
+viz(poly)
+viz!(other, color = :black, alpha = 0.2)
+viz!(boundary(clipped), color = :red, segmentsize = 3)
+Mke.current_figure()
+```
+
+## Weiler-Atherton
+
+```@docs
+WeilerAthertonClipping
+```
+
+```@example clipping
+# polygon to clip
+outer = Ring((8, 0), (4, 8), (2, 8), (-2, 0), (0, 0), (1, 2), (5, 2), (6, 0))
+inner = Ring((4, 4), (2, 4), (3, 6))
+poly = PolyArea([outer, inner])
+
+# clipping geometry can be a polygon with holes
+other = poly |> Rotate(π) |> Translate(8, 8)
+
+# clipped polygon
+clipped = clip(poly, other, WeilerAthertonClipping())
+
 viz(poly)
 viz!(other, color = :black, alpha = 0.2)
 viz!(boundary(clipped), color = :red, segmentsize = 3)
diff --git a/src/Meshes.jl b/src/Meshes.jl
index 58965ef93..9137fff34 100644
--- a/src/Meshes.jl
+++ b/src/Meshes.jl
@@ -440,6 +440,7 @@ export
   # clipping
   ClippingMethod,
   SutherlandHodgmanClipping,
+  WeilerAthertonClipping,
   clip,
 
   # intersections
diff --git a/src/clipping.jl b/src/clipping.jl
index aa9e1560d..0aab27182 100644
--- a/src/clipping.jl
+++ b/src/clipping.jl
@@ -21,3 +21,4 @@ function clip end
 # ----------------
 
 include("clipping/sutherlandhodgman.jl")
+include("clipping/weileratherton.jl")
diff --git a/src/clipping/weileratherton.jl b/src/clipping/weileratherton.jl
new file mode 100644
index 000000000..f9552f00d
--- /dev/null
+++ b/src/clipping/weileratherton.jl
@@ -0,0 +1,432 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+"""
+    WeilerAthertonClipping()
+
+The Wieler-Atherton algorithm for clipping polygons.
+
+## References
+
+* Weiler, K., & Atherton, P. 1977. [Hidden surface removal using polygon area sorting]
+  (https://dl.acm.org/doi/pdf/10.1145/563858.563896)
+"""
+struct WeilerAthertonClipping <: ClippingMethod end
+
+# VertexType has three different types of vertices used for deciding how to navigate the data
+# structure when collecting the polygon rings after clipping.
+abstract type VertexType end
+abstract type Normal <: VertexType end
+abstract type Entering <: VertexType end
+abstract type Exiting <: VertexType end
+
+# Data structure for clipping the polygons. Fields left and right are used depending on the
+# VertexType. If Normal, left points to the following vertex of the original ring, and right
+# to the next  intersection vertex on the edge between point and left.point, or the following
+# original ring vertex if no intersections on the edge. For Entering and Exiting types, left
+# poitns to the following vertex on the clipping ring and right to the following vertex on one
+# of the clipped rings.
+mutable struct RingVertex{VT<:VertexType,M<:Manifold,C<:CRS}
+  point::Point{M,C}
+  left::RingVertex{VT,M,C} where {VT<:VertexType}
+  right::RingVertex{VT,M,C} where {VT<:VertexType}
+
+  function RingVertex{VT,M,C}(point::Point{M,C}) where {VT<:VertexType,M<:Manifold,C<:CRS}
+    v = new(point)
+    v.left = v
+    v.right = v
+    return v
+  end
+end
+
+RingVertex{VT}(point::Point{M,C}) where {VT<:VertexType,M<:Manifold,C<:CRS} =
+  RingVertex{VT,manifold(point),typeof(point.coords)}(point)
+
+isnormal(::RingVertex{Normal}) = true
+isnormal(::RingVertex) = false
+
+function appendvertices!(v1::RingVertex{Normal}, v2::RingVertex{Normal})
+  v2.left = v1.left
+  v1.left = v2
+  v2.right = v1.right
+  v1.right = v2
+end
+
+function clip(poly::Polygon, ring::Ring, ::WeilerAthertonClipping)
+  polyrings = rings(poly)
+
+  # Convert the subject polygon rings and the clipping ring to the RingVertex data structure.
+  clippedrings = [gettraversalring(r) for r in polyrings]
+  startclipping = gettraversalring(ring)
+
+  # For keeping track of intersected rings, as the non-intersected ones need additional
+  # processing at the end.
+  intersected = zeros(Bool, length(clippedrings))
+
+  # For collecting all entering vertices, as they are used as starting points for collection
+  # of the output rings.
+  entering = RingVertex{Entering}[]
+
+  clipping = startclipping
+  while true
+    # Three consecutive clipping vertices are used to properly identify intersection vertex type
+    # for corner cases, so the clipping segment is constructed with the two following vertices
+    # after the current one.
+    clippingsegment = Segment(clipping.left.point, clipping.left.left.point)
+
+    for (k, startclipped) in enumerate(clippedrings)
+      clipped = startclipped
+      while true
+        # Like for the clipping, the clipped also uses three consecutive vertices.
+        clippedsegment = Segment(clipped.left.point, clipped.left.left.point)
+
+        vertextype = nothing
+        I = intersection(clippingsegment, clippedsegment)
+
+        # First try to handle Crossing, EdgeTouching and CornerTouching intersections as they might
+        # add only a single intersection.
+
+        if type(I) == Crossing
+          point = get(I)
+          cl = Line(clipping.left.point, clipping.left.left.point)
+          vertextype = sideof(clipped.left.left.point, cl) == LEFT ? Entering : Exiting
+        elseif type(I) == EdgeTouching
+          point = get(I)
+          if point ≈ clipped.left.point
+            # When intersection is at the shared vertex of two edges of the clipped ring,
+            # then split the interscting edge of the clipping ring at the intersection point.
+            vertextype = decidedirection(
+              Segment(clipped.point, clipped.left.point),
+              clippedsegment,
+              Segment(clipping.left.point, point),
+              Segment(point, clipping.left.left.point)
+            )
+          elseif point ≈ clipping.left.point
+            # When intersection is at the shared vertex of two edges of the clipping ring,
+            # then split the interscting edge of the clipped ring at the intersection point.
+            vertextype = decidedirection(
+              Segment(clipped.left.point, point),
+              Segment(point, clipped.left.left.point),
+              Segment(clipping.point, clipping.left.point),
+              clippingsegment
+            )
+          end
+        elseif type(I) == CornerTouching
+          # When intersection is at the shared vertices of both the clipping and the clipped rings.
+          point = get(I)
+          if (point ≈ clipped.left.point) && (point ≈ clipping.left.point)
+            # Only applies if the intersection coincides with the middles of the currently observed
+            # vertices for both clipping and clipped rings.
+            vertextype = decidedirection(
+              Segment(clipped.point, point),
+              Segment(point, clipped.left.left.point),
+              Segment(clipping.point, point),
+              Segment(point, clipping.left.left.point)
+            )
+          end
+        end
+        if !isnothing(vertextype)
+          insertintersections!(clipping.left, clipped.left, point, vertextype, entering)
+          if vertextype == Entering
+            intersected[k] = true
+          end
+        end
+
+        # Overlapping intersections might add up to two intersections, so handle separately.
+
+        if type(I) == Overlapping
+          for point in extrema(get(I))
+            # For both ends of the intersecting segment, check if it coincides with the middle of
+            # the obeserved vertices for clipped and clipping rings. If it does, attempt adding a
+            # point.
+
+            if point ≈ clipped.left.point
+              clippingprev = Segment(clipping.left.point, point)
+              if measure(clippingprev) ≈ 0.0 * u"m"
+                clippingprev = Segment(clipping.point, point)
+              end
+
+              vertextype = decidedirection(
+                Segment(clipped.point, point),
+                Segment(point, clipped.left.left.point),
+                clippingprev,
+                Segment(point, clipping.left.left.point)
+              )
+              if !isnothing(vertextype)
+                insertintersections!(clipping.left, clipped.left, point, vertextype, entering)
+                if vertextype == Entering
+                  intersected[k] = true
+                end
+              end
+            end
+
+            if point ≈ clipping.left.point
+              clippedprev = Segment(clipped.left.point, point)
+              if measure(clippedprev) ≈ 0.0 * u"m"
+                clippedprev = Segment(clipped.point, point)
+              end
+
+              vertextype = decidedirection(
+                clippedprev,
+                Segment(point, clipped.left.left.point),
+                Segment(clipping.point, point),
+                Segment(point, clipping.left.left.point)
+              )
+              if !isnothing(vertextype)
+                insertintersections!(clipping.left, clipped.left, point, vertextype, entering)
+                if vertextype == Entering
+                  intersected[k] = true
+                end
+              end
+            end
+          end
+        end
+
+        clipped = clipped.left
+        if clipped.left == startclipped.left
+          break
+        end
+      end
+    end
+
+    clipping = clipping.left
+    if clipping.left == startclipping.left
+      break
+    end
+  end
+
+  collectedrings = collectclipped(entering)
+
+  # When no interesction have been registered with any of the rings, the clipping ring either
+  # encircles everything or is completely contained in the clipped polygon.
+  if all(isequal(false), intersected)
+    o = orientation(ring)
+    contained = [v ∈ PolyArea(ring) for v in vertices(polyrings[1])]
+    if (o == CCW && all(contained))
+      return poly
+    end
+    if (o == CW && !all(contained))
+      push!(collectedrings, polyrings[1])
+      intersected[1] = true
+      for polyring in polyrings[2:end]
+        if !all([v ∈ PolyArea(ring) for v in vertices(polyring)])
+          push!(collectedrings, polyring)
+        end
+      end
+      if any([v ∈ PolyArea(polyrings[1]) for v in vertices(ring)])
+        push!(collectedrings, ring)
+      end
+      return PolyArea(collectedrings...)
+    end
+    if all([v ∈ PolyArea(polyrings[1]) for v in vertices(ring)])
+      push!(collectedrings, ring)
+    end
+  end
+
+  # Handle all the non-intersected rings. Add them if they are contained in the clipping ring.
+  polys = PolyArea[]
+  for r in collectedrings
+    newpolyrings = [r]
+    valid = true
+    for k in eachindex(intersected)
+      if !intersected[k]
+        if orientation(ring) == CCW
+          # Check if majority of vertices are inside the clipping ring.
+          ins = count(isequal(true), [v ∈ PolyArea(ring) for v in vertices(polyrings[k])])
+          if ins > (length(vertices(polyrings[k])) - ins)
+            if orientation(polyrings[k]) == CCW
+              pushfirst!(newpolyrings, polyrings[k])
+            else
+              push!(newpolyrings, polyrings[k])
+            end
+            intersected[k] = true
+          elseif vertices(ring)[1] ∈ PolyArea(polyrings[k]) && k > 1
+            # If the ring is contained within one of the inner rings, invalidate it.
+            valid = false
+          end
+        else
+          # Check if majority of vertices are outside the clipping ring.
+          ins = count(isequal(true), [v ∈ PolyArea(ring) for v in vertices(polyrings[k])])
+          if ins < (length(vertices(polyrings[k])) - ins)
+            if orientation(polyrings[k]) == CCW
+              pushfirst!(newpolyrings, polyrings[k])
+            else
+              push!(newpolyrings, polyrings[k])
+            end
+            intersected[k] = true
+          end
+        end
+      end
+    end
+    if valid
+      push!(polys, PolyArea(newpolyrings))
+    end
+  end
+
+  n = length(polys)
+  out = n == 0 ? nothing : (n == 1 ? polys[1] : GeometrySet(polys))
+  return out
+end
+
+function clip(poly::Polygon, other::Geometry, method::WeilerAthertonClipping)
+  return _clip(poly, boundary(other), method)
+end
+
+function _clip(poly, multi::Multi, method::WeilerAthertonClipping)
+  for r in parent(multi)
+    poly = clip(poly, r, method)
+    if isnothing(poly)
+      return nothing
+    end
+  end
+  return poly
+end
+
+_clip(poly, other, method) = clip(poly, other, method::WeilerAthertonClipping)
+
+function clip(dom::Domain, other::Geometry, method::WeilerAthertonClipping)
+  clipped = filter(!isnothing, [clip(geom, other, method) for geom in dom])
+  return isempty(clipped) ? nothing : (length(clipped) == 1 ? clipped[1] : GeometrySet(clipped))
+end
+
+# Inserts the intersection in the ring.
+function insertintersection!(head::RingVertex, intersection::RingVertex, side::Symbol)
+  tail = head.left
+  vertex = head
+
+  new = measure(Segment(head.point, intersection.point))
+  while true
+    os = isnormal(vertex) ? :right : side
+    current = measure(Segment(head.point, getfield(vertex, os).point))
+    if (new < current) || (getfield(vertex, os) == tail)
+      next = getfield(vertex, os)
+      if !(new ≈ measure(Segment(head.point, next.point)))
+        setfield!(intersection, side, next)
+        setfield!(vertex, os, intersection)
+      end
+      break
+    end
+    vertex = getfield(vertex, os)
+  end
+end
+
+# Inserts the intersection into both the clipping and the clipped rings.
+function insertintersections!(clipping, clipped, point, vtype, entering)
+  intersection = RingVertex{vtype}(point)
+  insertintersection!(clipping, intersection, :left)
+  insertintersection!(clipped, intersection, :right)
+
+  if vtype == Entering
+    push!(entering, intersection)
+  end
+end
+
+# Traversing and selecting following vertices for collecting of the rings after clipping depends
+# on the vertex type. Helper nextvertex follows the correct path.
+nextvertex(v::RingVertex{Normal}) = v.right
+nextvertex(v::RingVertex{Entering}) = v.right
+nextvertex(v::RingVertex{Exiting}) = v.left
+
+# Takes a list of entering vertices and returns all rings that contain those vertices.
+function collectclipped(entering::Vector{RingVertex{Entering}})
+  rings::Vector{Ring} = []
+  visited::Vector{RingVertex} = []
+  for i in eachindex(entering)
+    if entering[i] in visited
+      continue
+    end
+
+    ring::Vector{RingVertex} = []
+    vertex = entering[i]
+    while !(vertex in ring)
+      if vertex in visited
+        break
+      end
+
+      push!(ring, vertex)
+      push!(visited, vertex)
+      vertex = nextvertex(vertex)
+    end
+
+    # Remove duplicates.
+    newring::Vector{RingVertex} = [ring[1]]
+    for i in 2:length(ring)
+      if !(ring[i].point ≈ newring[end].point)
+        push!(newring, ring[i])
+      end
+    end
+    ring = newring
+
+    # Polygon might start several vertices after the first collected.
+    # This generally happens when there are overlapping edges that lead
+    # to several entering vertices without exiting in between. Then, the
+    # actual polygon is found by discarding the extra vertices before the
+    # proper loop.
+    k = findfirst(x -> ring[end].point == x.point, ring[1:(end - 1)])
+    if !isnothing(k)
+      ring = ring[(k + 1):end]
+    end
+
+    if length(ring) > 2
+      push!(rings, Ring([r.point for r in ring]...))
+    end
+  end
+  return rings
+end
+
+# Used to figure out the type of the vertex to add when intersection is other than the crossing
+# type. For input it takes four segments which all overlap on a single central vertex.
+# Function measures the angles formed between the segments and checks wether the two clipped segments
+# start and end in the same region or different regions, as separated by the clipping segments.
+function decidedirection(clipped₁, clipped₂, clipping₁, clipping₂)
+  # The input segments should all share one common vertex.
+  # Form vectors from the common vertex to other vertices.
+  a = minimum(clipped₁) - maximum(clipped₁)
+  b = maximum(clipped₂) - minimum(clipped₂)
+  c = minimum(clipping₁) - maximum(clipping₁)
+  d = maximum(clipping₂) - minimum(clipping₂)
+
+  if any(norm.([a, b, c, d]) .≈ 0.0 * u"m")
+    # An zero length input segment found, no need to calculate.
+    return nothing
+  end
+
+  β = mod(∠(a, b), 2pi * u"rad")
+  γ = mod(∠(a, c), 2pi * u"rad")
+  δ = mod(∠(a, d), 2pi * u"rad")
+
+  if isapprox(γ, 0.0, atol=atol(0.0)) || isapprox(γ, 2pi, atol=atol(0.0))
+    if δ < β < 2pi
+      return Entering
+    else
+      return Exiting
+    end
+  elseif isapprox(δ, 0.0, atol=atol(0.0)) || isapprox(δ, 2pi, atol=atol(0.0))
+    if γ < β < 2pi
+      return Exiting
+    else
+      return Entering
+    end
+  end
+  if γ < δ && (γ < β || isapprox(β, γ, atol=atol(0.0))) && (isapprox(β, δ, atol=atol(0.0)) || β < δ)
+    return Exiting
+  elseif δ < γ && (δ < β || isapprox(β, δ, atol=atol(0.0))) && (isapprox(β, γ, atol=atol(0.0)) || β < γ)
+    return Entering
+  end
+  return nothing
+end
+
+# Converts a regular Meshes.Ring into a ring formed with RingVertex data type.
+function gettraversalring(ring::Ring)
+  vs = vertices(ring)
+  start = RingVertex{Normal}(vs[1])
+
+  for v in vs[2:end]
+    new = RingVertex{Normal}(v)
+    appendvertices!(start, new)
+    start = new
+  end
+
+  return start.left
+end
diff --git a/test/clipping.jl b/test/clipping.jl
index 117f8a8ce..dbf838e76 100644
--- a/test/clipping.jl
+++ b/test/clipping.jl
@@ -74,3 +74,176 @@
   @test all(vertices(crings[1]) .≈ [cart(6, 4), cart(6, 5), cart(1, 2.5), cart(1, 1), cart(5, 2)])
   @test all(vertices(crings[2]) .≈ [cart(3.0, 3.0), cart(3.0, 3.5), cart(10 / 3, 11 / 3), cart(4.0, 3.0)])
 end
+
+@testitem "WeilerAtherton" setup = [Setup] begin
+  # triangle
+  poly = Triangle(cart(6, 2), cart(3, 5), cart(0, 2))
+  other = Quadrangle(cart(5, 0), cart(5, 4), cart(0, 4), cart(0, 0))
+  clipped = clip(poly, other, WeilerAthertonClipping())
+  @test issimple(clipped)
+  @test all(vertices(clipped) .≈ [cart(2, 4), cart(0, 2), cart(5, 2), cart(5, 3), cart(4, 4)])
+
+  # octagon
+  poly = Octagon(cart(8, -2), cart(8, 5), cart(2, 5), cart(4, 3), cart(6, 3), cart(4, 1), cart(2, 1), cart(2, -2))
+  other = Quadrangle(cart(5, 0), cart(5, 4), cart(0, 4), cart(0, 0))
+  clipped = clip(poly, other, WeilerAthertonClipping())
+  @test length(clipped) == 2
+  @test issimple(clipped[1])
+  @test issimple(clipped[2])
+  out = vertices.(clipped)
+  @test all(out[1] .≈ [cart(3, 4), cart(4, 3), cart(5, 3), cart(5, 4)])
+  @test all(out[2] .≈ [cart(5, 2), cart(4, 1), cart(2, 1), cart(2, 0), cart(5, 0)])
+
+  # inside
+  poly = Quadrangle(cart(1, 0), cart(1, 1), cart(0, 1), cart(0, 0))
+  other = Quadrangle(cart(5, 0), cart(5, 4), cart(0, 4), cart(0, 0))
+  clipped = clip(poly, other, WeilerAthertonClipping())
+  @test issimple(clipped)
+  @test clipped ≗ poly
+
+  # outside
+  poly = Quadrangle(cart(7, 6), cart(7, 7), cart(6, 7), cart(6, 6))
+  other = Quadrangle(cart(5, 0), cart(5, 4), cart(0, 4), cart(0, 0))
+  clipped = clip(poly, other, WeilerAthertonClipping())
+  @test isnothing(clipped)
+
+  # surrounded
+  poly = Hexagon(cart(0, 2), cart(-2, 2), cart(-2, 0), cart(0, -2), cart(2, -2), cart(2, 0))
+  other = Hexagon(cart(1, 0), cart(0, 1), cart(-1, 1), cart(-1, 0), cart(0, -1), cart(1, -1))
+  clipped = clip(poly, other, WeilerAthertonClipping())
+  @test issimple(clipped)
+  @test clipped ≗ other
+
+  # PolyArea with box
+  outer = Ring(cart(8, 0), cart(4, 8), cart(2, 8), cart(-2, 0), cart(0, 0), cart(1, 2), cart(5, 2), cart(6, 0))
+  inner = Ring(cart(4, 4), cart(2, 4), cart(3, 6))
+  poly = PolyArea([outer, inner])
+  other = Box(cart(0, 1), cart(3, 7))
+  clipped = clip(poly, other, WeilerAthertonClipping())
+  crings = rings(clipped)
+  @test issimple(clipped)
+  @test all(
+    vertices(crings[1]) .≈ [
+      cart(3.0, 4.0),
+      cart(2.0, 4.0),
+      cart(3.0, 6.0),
+      cart(3.0, 7.0),
+      cart(1.5, 7.0),
+      cart(0.0, 4.0),
+      cart(0.0, 1.0),
+      cart(0.5, 1.0),
+      cart(1.0, 2.0),
+      cart(3.0, 2.0)
+    ]
+  )
+
+  # PolyArea with outer ring outside and inner ring inside
+  outer = Ring(cart(8, 0), cart(2, 6), cart(-4, 0))
+  inner = Ring(cart(1, 3), cart(3, 3), cart(3, 1), cart(1, 1))
+  poly = PolyArea([outer, inner])
+  other = Quadrangle(cart(4, 4), cart(0, 4), cart(0, 0), cart(4, 0))
+  clipped = clip(poly, other, WeilerAthertonClipping())
+  @test !issimple(clipped)
+  crings = rings(clipped)
+  @test all(vertices(crings[1]) .≈ [cart(4, 0), cart(4, 4), cart(0, 4), cart(0, 0)])
+  @test crings[2] ≗ inner
+
+  # PolyArea with one inner ring inside `other` and another inner ring outside `other`
+  outer = Ring(cart(6, 4), cart(6, 7), cart(1, 6), cart(1, 1), cart(5, 2))
+  inner₁ = Ring(cart(3, 3), cart(3, 4), cart(4, 3))
+  inner₂ = Ring(cart(2, 5), cart(2, 6), cart(3, 5))
+  poly = PolyArea([outer, inner₁, inner₂])
+  other = PolyArea(Ring(cart(6, 1), cart(7, 2), cart(6, 5), cart(0, 2), cart(1, 1)))
+  clipped = clip(poly, other, WeilerAthertonClipping())
+  crings = rings(clipped)
+  @test issimple(clipped)
+  @test length(crings) == 1
+  @test all(
+    vertices(crings[1]) .≈ [
+      cart(1, 2.5),
+      cart(1, 1),
+      cart(5, 2),
+      cart(6, 4),
+      cart(6, 5),
+      cart(10 / 3, 11 / 3),
+      cart(4.0, 3.0),
+      cart(3.0, 3.0),
+      cart(3.0, 3.5)
+    ]
+  )
+
+  # Two polygons with holes.
+  outer = Ring((8, 0), (4, 8), (2, 8), (-2, 0), (0, 0), (1, 2), (5, 2), (6, 0))
+  inner = Ring((4, 4), (2, 4), (3, 6))
+  poly = PolyArea([outer, inner])
+  other = poly |> Rotate(pi) |> Translate(8.0, 8.0)
+  clipped = clip(poly, other, WeilerAthertonClipping())
+  crings = rings(clipped)
+  @test all(
+    vertices(crings[1]) .≈
+    [cart(7, 2), cart(5, 6), cart(3, 6), cart(2, 8), cart(1, 6), cart(3, 2), cart(5, 2), cart(6, 0)]
+  )
+  @test crings[2] ≗ inner
+  @test all(vertices(crings[3]) .≈ [cart(4, 4), cart(6, 4), cart(5, 2)])
+
+  # Tolerances are not properly retrieved for Float32 types, so need to pass them explicitly.
+  atol_el = coords(cart(0.0)).x
+
+  # Overlapping clipping polygon with the edges of the hole of the subject polygon.
+  triangle = Triangle((0.0, 0.0), (1.0, 0.0), (0.5, 1.0))
+  rectangle = Quadrangle((0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0))
+  poly = PolyArea(rings(rectangle)[1], reverse(rings(triangle |> Scale(0.5) |> Translate((0.3, 0.3)))[1]))
+  other = triangle |> Scale(0.9) |> Translate(0.1, -0.1)
+  clipped = clip(poly, other, WeilerAthertonClipping())
+  crings = rings(clipped)
+  @test length(crings) == 1
+  @test all(
+    isapprox.(
+      vertices(crings[1]),
+      [cart(0.8, 0.3), cart(0.3, 0.3), cart(0.15, 0.0), cart(0.95, 0.0)],
+      atol=atol(atol_el)
+    )
+  )
+
+  # Intersection only at a corner of two polygons with holes.
+  other = poly |> Rotate(pi)
+  clipped = clip(poly, other, WeilerAthertonClipping())
+  @test isnothing(clipped)
+
+  # Proper intersection of two polygons with holes.
+  other = poly |> Rotate(pi) |> Translate(1.0, 1.0)
+  clipped = clip(poly, other, WeilerAthertonClipping())
+  crings = rings(clipped)
+  @test length(crings) == 2
+  @test crings[1] ≈ rings(rectangle)[1]
+  @test all(
+    isapprox.(
+      vertices(crings[2]),
+      [
+        cart(0.4, 0.3),
+        cart(0.3, 0.3),
+        cart(0.35, 0.4),
+        cart(0.2, 0.7),
+        cart(0.5, 0.7),
+        cart(0.55, 0.8),
+        cart(0.6, 0.7),
+        cart(0.7, 0.7),
+        cart(0.65, 0.6),
+        cart(0.8, 0.3),
+        cart(0.5, 0.3),
+        cart(0.45, 0.2)
+      ],
+      atol=atol(atol_el)
+    )
+  )
+  # Clipping a GeometrySet.
+  poly = Quadrangle(cart(0.0, 0.0), cart(1.0, 0.0), cart(1.0, 1.0), cart(0.0, 1.0))
+  polyset = GeometrySet([poly, poly |> Translate(2.0, 0.0)])
+  other = Quadrangle(cart(0.5, 0.25), cart(2.5, 0.25), cart(2.5, 0.75), cart(0.5, 0.75))
+  clipped = clip(polyset, other, WeilerAthertonClipping())
+  @test clipped isa GeometrySet
+  crings = rings.(clipped)
+  @test length(crings) == 2
+  @test all(vertices(crings[1][1]) .≈ [cart(1.0, 0.25), cart(1.0, 0.75), cart(0.5, 0.75), cart(0.5, 0.25)])
+  @test all(vertices(crings[2][1]) .≈ [cart(2.0, 0.75), cart(2.0, 0.25), cart(2.5, 0.25), cart(2.5, 0.75)])
+end

From 3d0d0287f26f96460b98a309889bc38622faac3f Mon Sep 17 00:00:00 2001
From: Elias Carvalho <eliascarvdev@gmail.com>
Date: Wed, 4 Sep 2024 15:27:36 -0300
Subject: [PATCH 2/6] Code style updates

---
 src/clipping/weileratherton.jl | 101 ++++++++++++++++-----------------
 1 file changed, 50 insertions(+), 51 deletions(-)

diff --git a/src/clipping/weileratherton.jl b/src/clipping/weileratherton.jl
index f9552f00d..476ca860d 100644
--- a/src/clipping/weileratherton.jl
+++ b/src/clipping/weileratherton.jl
@@ -29,28 +29,27 @@ abstract type Exiting <: VertexType end
 # of the clipped rings.
 mutable struct RingVertex{VT<:VertexType,M<:Manifold,C<:CRS}
   point::Point{M,C}
-  left::RingVertex{VT,M,C} where {VT<:VertexType}
-  right::RingVertex{VT,M,C} where {VT<:VertexType}
+  left::RingVertex{<:VertexType,M,C}
+  right::RingVertex{<:VertexType,M,C}
 
-  function RingVertex{VT,M,C}(point::Point{M,C}) where {VT<:VertexType,M<:Manifold,C<:CRS}
+  function RingVertex{VT,M,C}(point) where {VT<:VertexType,M<:Manifold,C<:CRS}
     v = new(point)
     v.left = v
     v.right = v
-    return v
+    v
   end
 end
 
-RingVertex{VT}(point::Point{M,C}) where {VT<:VertexType,M<:Manifold,C<:CRS} =
-  RingVertex{VT,manifold(point),typeof(point.coords)}(point)
+RingVertex{VT}(point::Point{M,C}) where {VT<:VertexType,M<:Manifold,C<:CRS} = RingVertex{VT,M,C}(point)
 
 isnormal(::RingVertex{Normal}) = true
 isnormal(::RingVertex) = false
 
-function appendvertices!(v1::RingVertex{Normal}, v2::RingVertex{Normal})
-  v2.left = v1.left
-  v1.left = v2
-  v2.right = v1.right
-  v1.right = v2
+function appendvertices!(v₁::RingVertex{Normal}, v₂::RingVertex{Normal})
+  v₂.left = v₁.left
+  v₁.left = v₂
+  v₂.right = v₁.right
+  v₁.right = v₂
 end
 
 function clip(poly::Polygon, ring::Ring, ::WeilerAthertonClipping)
@@ -143,7 +142,7 @@ function clip(poly::Polygon, ring::Ring, ::WeilerAthertonClipping)
 
             if point ≈ clipped.left.point
               clippingprev = Segment(clipping.left.point, point)
-              if measure(clippingprev) ≈ 0.0 * u"m"
+              if measure(clippingprev) ≈ 0.0u"m"
                 clippingprev = Segment(clipping.point, point)
               end
 
@@ -163,7 +162,7 @@ function clip(poly::Polygon, ring::Ring, ::WeilerAthertonClipping)
 
             if point ≈ clipping.left.point
               clippedprev = Segment(clipped.left.point, point)
-              if measure(clippedprev) ≈ 0.0 * u"m"
+              if measure(clippedprev) ≈ 0.0u"m"
                 clippedprev = Segment(clipped.point, point)
               end
 
@@ -200,26 +199,26 @@ function clip(poly::Polygon, ring::Ring, ::WeilerAthertonClipping)
 
   # When no interesction have been registered with any of the rings, the clipping ring either
   # encircles everything or is completely contained in the clipped polygon.
-  if all(isequal(false), intersected)
+  if !any(intersected)
     o = orientation(ring)
-    contained = [v ∈ PolyArea(ring) for v in vertices(polyrings[1])]
-    if (o == CCW && all(contained))
+    allcontained = all(v ∈ PolyArea(ring) for v in vertices(polyrings[1]))
+    if (o == CCW && allcontained)
       return poly
     end
-    if (o == CW && !all(contained))
+    if (o == CW && !allcontained)
       push!(collectedrings, polyrings[1])
       intersected[1] = true
       for polyring in polyrings[2:end]
-        if !all([v ∈ PolyArea(ring) for v in vertices(polyring)])
+        if !all(v ∈ PolyArea(ring) for v in vertices(polyring))
           push!(collectedrings, polyring)
         end
       end
-      if any([v ∈ PolyArea(polyrings[1]) for v in vertices(ring)])
+      if any(v ∈ PolyArea(polyrings[1]) for v in vertices(ring))
         push!(collectedrings, ring)
       end
       return PolyArea(collectedrings...)
     end
-    if all([v ∈ PolyArea(polyrings[1]) for v in vertices(ring)])
+    if all(v ∈ PolyArea(polyrings[1]) for v in vertices(ring))
       push!(collectedrings, ring)
     end
   end
@@ -231,9 +230,9 @@ function clip(poly::Polygon, ring::Ring, ::WeilerAthertonClipping)
     valid = true
     for k in eachindex(intersected)
       if !intersected[k]
+        # Check if majority of vertices are inside the clipping ring.
+        ins = count(v ∈ PolyArea(ring) for v in vertices(polyrings[k]))
         if orientation(ring) == CCW
-          # Check if majority of vertices are inside the clipping ring.
-          ins = count(isequal(true), [v ∈ PolyArea(ring) for v in vertices(polyrings[k])])
           if ins > (length(vertices(polyrings[k])) - ins)
             if orientation(polyrings[k]) == CCW
               pushfirst!(newpolyrings, polyrings[k])
@@ -246,8 +245,6 @@ function clip(poly::Polygon, ring::Ring, ::WeilerAthertonClipping)
             valid = false
           end
         else
-          # Check if majority of vertices are outside the clipping ring.
-          ins = count(isequal(true), [v ∈ PolyArea(ring) for v in vertices(polyrings[k])])
           if ins < (length(vertices(polyrings[k])) - ins)
             if orientation(polyrings[k]) == CCW
               pushfirst!(newpolyrings, polyrings[k])
@@ -265,29 +262,24 @@ function clip(poly::Polygon, ring::Ring, ::WeilerAthertonClipping)
   end
 
   n = length(polys)
-  out = n == 0 ? nothing : (n == 1 ? polys[1] : GeometrySet(polys))
-  return out
+  n == 0 ? nothing : (n == 1 ? polys[1] : GeometrySet(polys))
 end
 
-function clip(poly::Polygon, other::Geometry, method::WeilerAthertonClipping)
-  return _clip(poly, boundary(other), method)
-end
+clip(poly::Polygon, other::Geometry, method::WeilerAthertonClipping) = clip(poly, boundary(other), method)
 
-function _clip(poly, multi::Multi, method::WeilerAthertonClipping)
+function clip(poly::Polygon, multi::Multi, method::WeilerAthertonClipping)
   for r in parent(multi)
     poly = clip(poly, r, method)
     if isnothing(poly)
       return nothing
     end
   end
-  return poly
+  poly
 end
 
-_clip(poly, other, method) = clip(poly, other, method::WeilerAthertonClipping)
-
 function clip(dom::Domain, other::Geometry, method::WeilerAthertonClipping)
   clipped = filter(!isnothing, [clip(geom, other, method) for geom in dom])
-  return isempty(clipped) ? nothing : (length(clipped) == 1 ? clipped[1] : GeometrySet(clipped))
+  isempty(clipped) ? nothing : (length(clipped) == 1 ? clipped[1] : GeometrySet(clipped))
 end
 
 # Inserts the intersection in the ring.
@@ -330,14 +322,14 @@ nextvertex(v::RingVertex{Exiting}) = v.left
 
 # Takes a list of entering vertices and returns all rings that contain those vertices.
 function collectclipped(entering::Vector{RingVertex{Entering}})
-  rings::Vector{Ring} = []
-  visited::Vector{RingVertex} = []
+  rings = Ring[]
+  visited = RingVertex[]
   for i in eachindex(entering)
     if entering[i] in visited
       continue
     end
 
-    ring::Vector{RingVertex} = []
+    ring = RingVertex[]
     vertex = entering[i]
     while !(vertex in ring)
       if vertex in visited
@@ -350,7 +342,7 @@ function collectclipped(entering::Vector{RingVertex{Entering}})
     end
 
     # Remove duplicates.
-    newring::Vector{RingVertex} = [ring[1]]
+    newring = RingVertex[ring[1]]
     for i in 2:length(ring)
       if !(ring[i].point ≈ newring[end].point)
         push!(newring, ring[i])
@@ -369,10 +361,10 @@ function collectclipped(entering::Vector{RingVertex{Entering}})
     end
 
     if length(ring) > 2
-      push!(rings, Ring([r.point for r in ring]...))
+      push!(rings, Ring([r.point for r in ring]))
     end
   end
-  return rings
+  rings
 end
 
 # Used to figure out the type of the vertex to add when intersection is other than the crossing
@@ -380,6 +372,8 @@ end
 # Function measures the angles formed between the segments and checks wether the two clipped segments
 # start and end in the same region or different regions, as separated by the clipping segments.
 function decidedirection(clipped₁, clipped₂, clipping₁, clipping₂)
+  T = numtype(lentype(clipped₁))
+
   # The input segments should all share one common vertex.
   # Form vectors from the common vertex to other vertices.
   a = minimum(clipped₁) - maximum(clipped₁)
@@ -387,34 +381,39 @@ function decidedirection(clipped₁, clipped₂, clipping₁, clipping₂)
   c = minimum(clipping₁) - maximum(clipping₁)
   d = maximum(clipping₂) - minimum(clipping₂)
 
-  if any(norm.([a, b, c, d]) .≈ 0.0 * u"m")
+  if any(norm(v) ≈ 0.0u"m" for v in (a, b, c, d))
     # An zero length input segment found, no need to calculate.
     return nothing
   end
 
-  β = mod(∠(a, b), 2pi * u"rad")
-  γ = mod(∠(a, c), 2pi * u"rad")
-  δ = mod(∠(a, d), 2pi * u"rad")
+  tol = atol(0.0) * u"rad"
+  twoπ = 2 * T(π) * u"rad"
 
-  if isapprox(γ, 0.0, atol=atol(0.0)) || isapprox(γ, 2pi, atol=atol(0.0))
-    if δ < β < 2pi
+  β = mod(∠(a, b), twoπ)
+  γ = mod(∠(a, c), twoπ)
+  δ = mod(∠(a, d), twoπ)
+
+  if isapprox(γ, zero(γ), atol=tol) || isapprox(γ, twoπ, atol=tol)
+    if δ < β < twoπ
       return Entering
     else
       return Exiting
     end
-  elseif isapprox(δ, 0.0, atol=atol(0.0)) || isapprox(δ, 2pi, atol=atol(0.0))
-    if γ < β < 2pi
+  elseif isapprox(δ, zero(δ), atol=tol) || isapprox(δ, twoπ, atol=tol)
+    if γ < β < twoπ
       return Exiting
     else
       return Entering
     end
   end
-  if γ < δ && (γ < β || isapprox(β, γ, atol=atol(0.0))) && (isapprox(β, δ, atol=atol(0.0)) || β < δ)
+
+  if γ < δ && (γ < β || isapprox(β, γ, atol=tol)) && (isapprox(β, δ, atol=tol) || β < δ)
     return Exiting
-  elseif δ < γ && (δ < β || isapprox(β, δ, atol=atol(0.0))) && (isapprox(β, γ, atol=atol(0.0)) || β < γ)
+  elseif δ < γ && (δ < β || isapprox(β, δ, atol=tol)) && (isapprox(β, γ, atol=tol) || β < γ)
     return Entering
   end
-  return nothing
+
+  nothing
 end
 
 # Converts a regular Meshes.Ring into a ring formed with RingVertex data type.

From 01ce2ebfde358790f442b62552916a2c4f29e8d3 Mon Sep 17 00:00:00 2001
From: Nemanja Jovanovic <nemanja.jovanovic@proton.me>
Date: Thu, 12 Sep 2024 11:14:23 +0300
Subject: [PATCH 3/6] Break up and refactor clip function and remove
 non-specific methods

---
 docs/src/algorithms/clipping.md |   4 +-
 src/clipping/weileratherton.jl  | 369 ++++++++++++++++----------------
 test/clipping.jl                |  66 +-----
 3 files changed, 185 insertions(+), 254 deletions(-)

diff --git a/docs/src/algorithms/clipping.md b/docs/src/algorithms/clipping.md
index 947126b58..2c0fd824d 100644
--- a/docs/src/algorithms/clipping.md
+++ b/docs/src/algorithms/clipping.md
@@ -46,8 +46,8 @@ outer = Ring((8, 0), (4, 8), (2, 8), (-2, 0), (0, 0), (1, 2), (5, 2), (6, 0))
 inner = Ring((4, 4), (2, 4), (3, 6))
 poly = PolyArea([outer, inner])
 
-# clipping geometry can be a polygon with holes
-other = poly |> Rotate(π) |> Translate(8, 8)
+# clipping geometry
+other = outer |> Rotate(π) |> Translate(8, 8)
 
 # clipped polygon
 clipped = clip(poly, other, WeilerAthertonClipping())
diff --git a/src/clipping/weileratherton.jl b/src/clipping/weileratherton.jl
index 476ca860d..5c11630d3 100644
--- a/src/clipping/weileratherton.jl
+++ b/src/clipping/weileratherton.jl
@@ -23,10 +23,11 @@ abstract type Exiting <: VertexType end
 
 # Data structure for clipping the polygons. Fields left and right are used depending on the
 # VertexType. If Normal, left points to the following vertex of the original ring, and right
-# to the next  intersection vertex on the edge between point and left.point, or the following
+# to the next intersection vertex on the edge between point and left.point, or the following
 # original ring vertex if no intersections on the edge. For Entering and Exiting types, left
 # poitns to the following vertex on the clipping ring and right to the following vertex on one
 # of the clipped rings.
+# TODO Either properly document the usage of left and right, or separate RingVertex into Entering, Exiting and Normal vertices.
 mutable struct RingVertex{VT<:VertexType,M<:Manifold,C<:CRS}
   point::Point{M,C}
   left::RingVertex{<:VertexType,M,C}
@@ -52,9 +53,21 @@ function appendvertices!(v₁::RingVertex{Normal}, v₂::RingVertex{Normal})
   v₁.right = v₂
 end
 
+# Traversing and selecting following vertices for collecting of the rings after clipping depends
+# on the vertex type. Helper nextvertex follows the correct path.
+nextvertex(v::RingVertex{Normal}) = v.right
+nextvertex(v::RingVertex{Entering}) = v.right
+nextvertex(v::RingVertex{Exiting}) = v.left
+
 function clip(poly::Polygon, ring::Ring, ::WeilerAthertonClipping)
   polyrings = rings(poly)
 
+  # If the polygon is contained in the ring, return the polygon right away.
+  allcontained = all(v ∈ PolyArea(ring) for v in vertices(polyrings[1]))
+  if (orientation(ring) == CCW && allcontained)
+    return poly
+  end
+
   # Convert the subject polygon rings and the clipping ring to the RingVertex data structure.
   clippedrings = [gettraversalring(r) for r in polyrings]
   startclipping = gettraversalring(ring)
@@ -80,206 +93,88 @@ function clip(poly::Polygon, ring::Ring, ::WeilerAthertonClipping)
         # Like for the clipping, the clipped also uses three consecutive vertices.
         clippedsegment = Segment(clipped.left.point, clipped.left.left.point)
 
-        vertextype = nothing
         I = intersection(clippingsegment, clippedsegment)
-
-        # First try to handle Crossing, EdgeTouching and CornerTouching intersections as they might
-        # add only a single intersection.
-
-        if type(I) == Crossing
-          point = get(I)
-          cl = Line(clipping.left.point, clipping.left.left.point)
-          vertextype = sideof(clipped.left.left.point, cl) == LEFT ? Entering : Exiting
-        elseif type(I) == EdgeTouching
-          point = get(I)
-          if point ≈ clipped.left.point
-            # When intersection is at the shared vertex of two edges of the clipped ring,
-            # then split the interscting edge of the clipping ring at the intersection point.
-            vertextype = decidedirection(
-              Segment(clipped.point, clipped.left.point),
-              clippedsegment,
-              Segment(clipping.left.point, point),
-              Segment(point, clipping.left.left.point)
-            )
-          elseif point ≈ clipping.left.point
-            # When intersection is at the shared vertex of two edges of the clipping ring,
-            # then split the interscting edge of the clipped ring at the intersection point.
-            vertextype = decidedirection(
-              Segment(clipped.left.point, point),
-              Segment(point, clipped.left.left.point),
-              Segment(clipping.point, clipping.left.point),
-              clippingsegment
-            )
-          end
-        elseif type(I) == CornerTouching
-          # When intersection is at the shared vertices of both the clipping and the clipped rings.
-          point = get(I)
-          if (point ≈ clipped.left.point) && (point ≈ clipping.left.point)
-            # Only applies if the intersection coincides with the middles of the currently observed
-            # vertices for both clipping and clipped rings.
-            vertextype = decidedirection(
-              Segment(clipped.point, point),
-              Segment(point, clipped.left.left.point),
-              Segment(clipping.point, point),
-              Segment(point, clipping.left.left.point)
-            )
-          end
-        end
-        if !isnothing(vertextype)
-          insertintersections!(clipping.left, clipped.left, point, vertextype, entering)
-          if vertextype == Entering
-            intersected[k] = true
-          end
-        end
-
-        # Overlapping intersections might add up to two intersections, so handle separately.
-
-        if type(I) == Overlapping
-          for point in extrema(get(I))
-            # For both ends of the intersecting segment, check if it coincides with the middle of
-            # the obeserved vertices for clipped and clipping rings. If it does, attempt adding a
-            # point.
-
-            if point ≈ clipped.left.point
-              clippingprev = Segment(clipping.left.point, point)
-              if measure(clippingprev) ≈ 0.0u"m"
-                clippingprev = Segment(clipping.point, point)
-              end
-
-              vertextype = decidedirection(
-                Segment(clipped.point, point),
-                Segment(point, clipped.left.left.point),
-                clippingprev,
-                Segment(point, clipping.left.left.point)
-              )
-              if !isnothing(vertextype)
-                insertintersections!(clipping.left, clipped.left, point, vertextype, entering)
-                if vertextype == Entering
-                  intersected[k] = true
-                end
-              end
-            end
-
-            if point ≈ clipping.left.point
-              clippedprev = Segment(clipped.left.point, point)
-              if measure(clippedprev) ≈ 0.0u"m"
-                clippedprev = Segment(clipped.point, point)
-              end
-
-              vertextype = decidedirection(
-                clippedprev,
-                Segment(point, clipped.left.left.point),
-                Segment(clipping.point, point),
-                Segment(point, clipping.left.left.point)
-              )
-              if !isnothing(vertextype)
-                insertintersections!(clipping.left, clipped.left, point, vertextype, entering)
-                if vertextype == Entering
-                  intersected[k] = true
-                end
-              end
-            end
-          end
-        end
+        vertex = vertexfromintersection(I, clipping, clipped)
+        success = insertintersections!(vertex, clipping, clipped, entering)
+        intersected[k] = intersected[k] || success
 
         clipped = clipped.left
-        if clipped.left == startclipped.left
-          break
-        end
+        clipped.left == startclipped.left && break
       end
     end
 
     clipping = clipping.left
-    if clipping.left == startclipping.left
-      break
-    end
+    clipping.left == startclipping.left && break
   end
 
-  collectedrings = collectclipped(entering)
-
-  # When no interesction have been registered with any of the rings, the clipping ring either
-  # encircles everything or is completely contained in the clipped polygon.
+  # Handle the case when no interections have been found.
   if !any(intersected)
-    o = orientation(ring)
-    allcontained = all(v ∈ PolyArea(ring) for v in vertices(polyrings[1]))
-    if (o == CCW && allcontained)
-      return poly
-    end
-    if (o == CW && !allcontained)
-      push!(collectedrings, polyrings[1])
-      intersected[1] = true
-      for polyring in polyrings[2:end]
-        if !all(v ∈ PolyArea(ring) for v in vertices(polyring))
-          push!(collectedrings, polyring)
-        end
-      end
-      if any(v ∈ PolyArea(polyrings[1]) for v in vertices(ring))
-        push!(collectedrings, ring)
-      end
-      return PolyArea(collectedrings...)
-    end
-    if all(v ∈ PolyArea(polyrings[1]) for v in vertices(ring))
-      push!(collectedrings, ring)
+    if orientation(ring) == CW
+      # For an inner clipping ring take all outside rings.
+      return PolyArea(collectoutsiderings(ring, polyrings)...)
+    else
+      # For an outer clipping ring add it to act as the outer ring.
+      collectedrings = all(v ∈ PolyArea(polyrings[1]) for v in vertices(ring)) ? [ring] : []
     end
+  else
+    # Collect rings formed from the intersected rings.
+    collectedrings = collectclipped(entering)
   end
 
-  # Handle all the non-intersected rings. Add them if they are contained in the clipping ring.
-  polys = PolyArea[]
-  for r in collectedrings
-    newpolyrings = [r]
-    valid = true
-    for k in eachindex(intersected)
-      if !intersected[k]
-        # Check if majority of vertices are inside the clipping ring.
-        ins = count(v ∈ PolyArea(ring) for v in vertices(polyrings[k]))
-        if orientation(ring) == CCW
-          if ins > (length(vertices(polyrings[k])) - ins)
-            if orientation(polyrings[k]) == CCW
-              pushfirst!(newpolyrings, polyrings[k])
-            else
-              push!(newpolyrings, polyrings[k])
-            end
-            intersected[k] = true
-          elseif vertices(ring)[1] ∈ PolyArea(polyrings[k]) && k > 1
-            # If the ring is contained within one of the inner rings, invalidate it.
-            valid = false
-          end
-        else
-          if ins < (length(vertices(polyrings[k])) - ins)
-            if orientation(polyrings[k]) == CCW
-              pushfirst!(newpolyrings, polyrings[k])
-            else
-              push!(newpolyrings, polyrings[k])
-            end
-            intersected[k] = true
-          end
-        end
-      end
-    end
-    if valid
-      push!(polys, PolyArea(newpolyrings))
-    end
-  end
+  # Complete the collected rings by adding non-intersected rings
+  # if they are contained within the collected ones.
+  completedrings = [perhapsaddnonintersected!(r, polyrings, intersected) for r in collectedrings]
+
+  # Convert completed ring lists into PolyAreas.
+  polys = PolyArea.(filter(!isnothing, completedrings))
 
   n = length(polys)
-  n == 0 ? nothing : (n == 1 ? polys[1] : GeometrySet(polys))
+  n == 0 ? nothing : (n == 1 ? polys[1] : Multi(polys))
 end
 
 clip(poly::Polygon, other::Geometry, method::WeilerAthertonClipping) = clip(poly, boundary(other), method)
 
-function clip(poly::Polygon, multi::Multi, method::WeilerAthertonClipping)
-  for r in parent(multi)
-    poly = clip(poly, r, method)
-    if isnothing(poly)
-      return nothing
+function collectoutsiderings(ring, polyrings)
+  newpolyrings = Ring[]
+  for k in eachindex(polyrings)
+    if !all(v ∈ PolyArea(ring) for v in vertices(polyrings[k]))
+      push!(newpolyrings, polyrings[k])
     end
   end
-  poly
+  if any(v ∈ PolyArea(polyrings[1]) for v in vertices(ring))
+    push!(newpolyrings, ring)
+  end
+  newpolyrings
 end
 
-function clip(dom::Domain, other::Geometry, method::WeilerAthertonClipping)
-  clipped = filter(!isnothing, [clip(geom, other, method) for geom in dom])
-  isempty(clipped) ? nothing : (length(clipped) == 1 ? clipped[1] : GeometrySet(clipped))
+function perhapsaddnonintersected!(ring, polyrings, intersected)
+  # Handle all the non-intersected rings. Add them if they are contained in the clipping ring.
+  newpolyrings = [ring]
+  for k in eachindex(intersected)
+    if !intersected[k]
+      ccw = orientation(ring) == CCW
+
+      # Discard if the processed ring is an outer ring but inside another inner ring.
+      if ccw && vertices(ring)[1] ∈ PolyArea(polyrings[k]) && k > 1
+        return nothing
+      end
+
+      vs = vertices(polyrings[k])
+      ins = count(v ∈ PolyArea(ring) for v in vs)
+      outs = length(vs) - ins
+      if ccw == (ins > outs)
+        if !ccw
+          # An outer ring should be first in the list.
+          pushfirst!(newpolyrings, polyrings[k])
+        else
+          # Inner rings should follow an outer in the list.
+          push!(newpolyrings, polyrings[k])
+        end
+        intersected[k] = true
+      end
+    end
+  end
+  newpolyrings
 end
 
 # Inserts the intersection in the ring.
@@ -304,21 +199,23 @@ function insertintersection!(head::RingVertex, intersection::RingVertex, side::S
 end
 
 # Inserts the intersection into both the clipping and the clipped rings.
-function insertintersections!(clipping, clipped, point, vtype, entering)
-  intersection = RingVertex{vtype}(point)
-  insertintersection!(clipping, intersection, :left)
-  insertintersection!(clipped, intersection, :right)
-
-  if vtype == Entering
-    push!(entering, intersection)
+function insertintersections!(vertex::Tuple, clipping, clipped, entering)
+  (vtype, point) = vertex
+  if !isnothing(vtype)
+    intersection = RingVertex{vtype}(point)
+    insertintersection!(clipping.left, intersection, :left)
+    insertintersection!(clipped.left, intersection, :right)
+
+    if vtype == Entering
+      push!(entering, intersection)
+    end
+    return true
   end
+  false
 end
 
-# Traversing and selecting following vertices for collecting of the rings after clipping depends
-# on the vertex type. Helper nextvertex follows the correct path.
-nextvertex(v::RingVertex{Normal}) = v.right
-nextvertex(v::RingVertex{Entering}) = v.right
-nextvertex(v::RingVertex{Exiting}) = v.left
+insertintersections!(vertices::Array, clipping, clipped, entering) =
+  any(insertintersections!.(vertices, Ref(clipping), Ref(clipped), Ref(entering)))
 
 # Takes a list of entering vertices and returns all rings that contain those vertices.
 function collectclipped(entering::Vector{RingVertex{Entering}})
@@ -367,6 +264,98 @@ function collectclipped(entering::Vector{RingVertex{Entering}})
   rings
 end
 
+function vertexfromintersection(I, clipping, clipped)
+  type(I) == Crossing && return vertexfromcrossing(get(I), clipping, clipped)
+  type(I) == CornerTouching && return vertexfromcornertouching(get(I), clipping, clipped)
+  type(I) == EdgeTouching && return vertexfromedgetouching(get(I), clipping, clipped)
+  type(I) == Overlapping && return vertexfromoverlapping(get(I), clipping, clipped)
+  (nothing, nothing)
+end
+
+function vertexfromcrossing(point, clipping, clipped)
+  cl = Line(clipping.left.point, clipping.left.left.point)
+  vertextype = sideof(clipped.left.left.point, cl) == LEFT ? Entering : Exiting
+  (vertextype, point)
+end
+
+function vertexfromedgetouching(point, clipping, clipped)
+  vertextype = nothing
+  if point ≈ clipped.left.point
+    # When intersection is at the shared vertex of two edges of the clipped ring,
+    # then split the interscting edge of the clipping ring at the intersection point.
+    vertextype = decidedirection(
+      Segment(clipped.point, clipped.left.point),
+      Segment(clipped.left.point, clipped.left.left.point),
+      Segment(clipping.left.point, point),
+      Segment(point, clipping.left.left.point)
+    )
+  elseif point ≈ clipping.left.point
+    # When intersection is at the shared vertex of two edges of the clipping ring,
+    # then split the interscting edge of the clipped ring at the intersection point.
+    vertextype = decidedirection(
+      Segment(clipped.left.point, point),
+      Segment(point, clipped.left.left.point),
+      Segment(clipping.point, clipping.left.point),
+      Segment(clipping.left.point, clipping.left.left.point)
+    )
+  end
+  (vertextype, point)
+end
+
+function vertexfromcornertouching(point, clipping, clipped)
+  vertextype = nothing
+  # When intersection is at the shared vertices of both the clipping and the clipped rings.
+  if (point ≈ clipped.left.point) && (point ≈ clipping.left.point)
+    # Only applies if the intersection coincides with the middles of the currently observed
+    # vertices for both clipping and clipped rings.
+    vertextype = decidedirection(
+      Segment(clipped.point, point),
+      Segment(point, clipped.left.left.point),
+      Segment(clipping.point, point),
+      Segment(point, clipping.left.left.point)
+    )
+  end
+  (vertextype, point)
+end
+
+function vertexfromoverlapping(segment, clipping, clipped)
+  # For both ends of the intersecting segment, check if it coincides with the middle of
+  # the obeserved vertices for clipped and clipping rings. If it does, attempt adding a
+  # point.
+
+  ret = Tuple[]
+  for point in extrema(segment)
+    if point ≈ clipped.left.point
+      clippingprev = Segment(clipping.left.point, point)
+      if measure(clippingprev) ≈ 0.0u"m"
+        clippingprev = Segment(clipping.point, point)
+      end
+      vertextype = decidedirection(
+        Segment(clipped.point, point),
+        Segment(point, clipped.left.left.point),
+        clippingprev,
+        Segment(point, clipping.left.left.point)
+      )
+      push!(ret, (vertextype, point))
+    end
+
+    if point ≈ clipping.left.point
+      clippedprev = Segment(clipped.left.point, point)
+      if measure(clippedprev) ≈ 0.0u"m"
+        clippedprev = Segment(clipped.point, point)
+      end
+      vertextype = decidedirection(
+        clippedprev,
+        Segment(point, clipped.left.left.point),
+        Segment(clipping.point, point),
+        Segment(point, clipping.left.left.point)
+      )
+      push!(ret, (vertextype, point))
+    end
+  end
+  ret
+end
+
 # Used to figure out the type of the vertex to add when intersection is other than the crossing
 # type. For input it takes four segments which all overlap on a single central vertex.
 # Function measures the angles formed between the segments and checks wether the two clipped segments
diff --git a/test/clipping.jl b/test/clipping.jl
index dbf838e76..7206afb75 100644
--- a/test/clipping.jl
+++ b/test/clipping.jl
@@ -87,12 +87,10 @@ end
   poly = Octagon(cart(8, -2), cart(8, 5), cart(2, 5), cart(4, 3), cart(6, 3), cart(4, 1), cart(2, 1), cart(2, -2))
   other = Quadrangle(cart(5, 0), cart(5, 4), cart(0, 4), cart(0, 0))
   clipped = clip(poly, other, WeilerAthertonClipping())
-  @test length(clipped) == 2
-  @test issimple(clipped[1])
-  @test issimple(clipped[2])
-  out = vertices.(clipped)
-  @test all(out[1] .≈ [cart(3, 4), cart(4, 3), cart(5, 3), cart(5, 4)])
-  @test all(out[2] .≈ [cart(5, 2), cart(4, 1), cart(2, 1), cart(2, 0), cart(5, 0)])
+  @test all(
+    vertices(clipped) .≈
+    [cart(3, 4), cart(4, 3), cart(5, 3), cart(5, 4), cart(5, 2), cart(4, 1), cart(2, 1), cart(2, 0), cart(5, 0)]
+  )
 
   # inside
   poly = Quadrangle(cart(1, 0), cart(1, 1), cart(0, 1), cart(0, 0))
@@ -172,20 +170,6 @@ end
     ]
   )
 
-  # Two polygons with holes.
-  outer = Ring((8, 0), (4, 8), (2, 8), (-2, 0), (0, 0), (1, 2), (5, 2), (6, 0))
-  inner = Ring((4, 4), (2, 4), (3, 6))
-  poly = PolyArea([outer, inner])
-  other = poly |> Rotate(pi) |> Translate(8.0, 8.0)
-  clipped = clip(poly, other, WeilerAthertonClipping())
-  crings = rings(clipped)
-  @test all(
-    vertices(crings[1]) .≈
-    [cart(7, 2), cart(5, 6), cart(3, 6), cart(2, 8), cart(1, 6), cart(3, 2), cart(5, 2), cart(6, 0)]
-  )
-  @test crings[2] ≗ inner
-  @test all(vertices(crings[3]) .≈ [cart(4, 4), cart(6, 4), cart(5, 2)])
-
   # Tolerances are not properly retrieved for Float32 types, so need to pass them explicitly.
   atol_el = coords(cart(0.0)).x
 
@@ -204,46 +188,4 @@ end
       atol=atol(atol_el)
     )
   )
-
-  # Intersection only at a corner of two polygons with holes.
-  other = poly |> Rotate(pi)
-  clipped = clip(poly, other, WeilerAthertonClipping())
-  @test isnothing(clipped)
-
-  # Proper intersection of two polygons with holes.
-  other = poly |> Rotate(pi) |> Translate(1.0, 1.0)
-  clipped = clip(poly, other, WeilerAthertonClipping())
-  crings = rings(clipped)
-  @test length(crings) == 2
-  @test crings[1] ≈ rings(rectangle)[1]
-  @test all(
-    isapprox.(
-      vertices(crings[2]),
-      [
-        cart(0.4, 0.3),
-        cart(0.3, 0.3),
-        cart(0.35, 0.4),
-        cart(0.2, 0.7),
-        cart(0.5, 0.7),
-        cart(0.55, 0.8),
-        cart(0.6, 0.7),
-        cart(0.7, 0.7),
-        cart(0.65, 0.6),
-        cart(0.8, 0.3),
-        cart(0.5, 0.3),
-        cart(0.45, 0.2)
-      ],
-      atol=atol(atol_el)
-    )
-  )
-  # Clipping a GeometrySet.
-  poly = Quadrangle(cart(0.0, 0.0), cart(1.0, 0.0), cart(1.0, 1.0), cart(0.0, 1.0))
-  polyset = GeometrySet([poly, poly |> Translate(2.0, 0.0)])
-  other = Quadrangle(cart(0.5, 0.25), cart(2.5, 0.25), cart(2.5, 0.75), cart(0.5, 0.75))
-  clipped = clip(polyset, other, WeilerAthertonClipping())
-  @test clipped isa GeometrySet
-  crings = rings.(clipped)
-  @test length(crings) == 2
-  @test all(vertices(crings[1][1]) .≈ [cart(1.0, 0.25), cart(1.0, 0.75), cart(0.5, 0.75), cart(0.5, 0.25)])
-  @test all(vertices(crings[2][1]) .≈ [cart(2.0, 0.75), cart(2.0, 0.25), cart(2.5, 0.25), cart(2.5, 0.75)])
 end

From 4dc59ca1d077ec4212bd96249ceae83fe564d3ea Mon Sep 17 00:00:00 2001
From: Nemanja Jovanovic <nemanja.jovanovic@proton.me>
Date: Tue, 12 Nov 2024 19:08:06 +0200
Subject: [PATCH 4/6] Use helper functions for RingVertex

Improve comments and algorithm description
---
 src/clipping/weileratherton.jl | 171 ++++++++++++++++++++-------------
 1 file changed, 104 insertions(+), 67 deletions(-)

diff --git a/src/clipping/weileratherton.jl b/src/clipping/weileratherton.jl
index 5c11630d3..1cf04a22a 100644
--- a/src/clipping/weileratherton.jl
+++ b/src/clipping/weileratherton.jl
@@ -14,20 +14,20 @@ The Wieler-Atherton algorithm for clipping polygons.
 """
 struct WeilerAthertonClipping <: ClippingMethod end
 
-# VertexType has three different types of vertices used for deciding how to navigate the data
-# structure when collecting the polygon rings after clipping.
+# VertexType distinguishes three different types of vertices used for constructing the structure
+# used by the Weiler-Atherton clipping algorithm. Normal type represents regular vertices, the
+# connecting points of edges of a rings. Entering and Exiting types represent intersections of
+# two ring edges, where one belongs to a clipping and other belongs to a clipped ring. The
+# Entering and Exiting types describe, from the perspective of the edges of the clipping ring,
+# whether it enters or exits the interior of the clipped ring.
 abstract type VertexType end
 abstract type Normal <: VertexType end
 abstract type Entering <: VertexType end
 abstract type Exiting <: VertexType end
 
 # Data structure for clipping the polygons. Fields left and right are used depending on the
-# VertexType. If Normal, left points to the following vertex of the original ring, and right
-# to the next intersection vertex on the edge between point and left.point, or the following
-# original ring vertex if no intersections on the edge. For Entering and Exiting types, left
-# poitns to the following vertex on the clipping ring and right to the following vertex on one
-# of the clipped rings.
-# TODO Either properly document the usage of left and right, or separate RingVertex into Entering, Exiting and Normal vertices.
+# VertexType, as designated with the helper functions. The data structure forms a directed
+# graph with each element always pointing to two elements.
 mutable struct RingVertex{VT<:VertexType,M<:Manifold,C<:CRS}
   point::Point{M,C}
   left::RingVertex{<:VertexType,M,C}
@@ -53,11 +53,41 @@ function appendvertices!(v₁::RingVertex{Normal}, v₂::RingVertex{Normal})
   v₁.right = v₂
 end
 
-# Traversing and selecting following vertices for collecting of the rings after clipping depends
-# on the vertex type. Helper nextvertex follows the correct path.
+# Helper functions to designate the use of the left and right branches of the RingVertex.
+
+# Normal vertex uses the left branch for the next Normal vertex,
+# which is the end of the ring edge for which the current vertex is the edge start.
+nextnormal(v::RingVertex{Normal}) = v.left
+# Normal vertex uses the right branch for the next intersection on the ring edge if any.
+# Otherwise, it holds the next Normal vertex, same as for the left branch.
 nextvertex(v::RingVertex{Normal}) = v.right
-nextvertex(v::RingVertex{Entering}) = v.right
-nextvertex(v::RingVertex{Exiting}) = v.left
+# Next clipping vertex is a following vertex, normal or an intersection, on a clipping ring.
+# Next clipped vertex is a following vertex, normal or an intersection, on a clipped ring.
+# Normal vertices are either on a clipping or a clipped, but the implementation is same
+# and helper functions can be used interchangeably.
+getnextclippingvertex(v::RingVertex{Normal}) = v.right
+setnextclippingvertex!(v::RingVertex{Normal}, next::RingVertex) = v.right = next
+getnextclippedvertex(v::RingVertex{Normal}) = v.right
+setnextclippedvertex!(v::RingVertex{Normal}, next::RingVertex) = v.right = next
+
+# Both Entering and Exiting vertices are always intersections and use the left branch for
+# the following vertex on the clipping ring.
+getnextclippingvertex(v::RingVertex) = v.left
+setnextclippingvertex!(v::RingVertex, next::RingVertex) = v.left = next
+
+# Similarly, the right branch holds the following vertex on the clipped ring.
+getnextclippedvertex(v::RingVertex) = v.right
+setnextclippedvertex!(v::RingVertex, next::RingVertex) = v.right = next
+
+# Traversing and selecting correct following vertices from the directed graph of RingVertex
+# elements should switch between rings whenever an intersection is encountered. If the current
+# edge is entering the interior of the clipped ring, follow to the next vertex on the clipping
+# ring. In case the edge is exiting the interior of the clipped ring, then follow the next vertex
+# on the clipped ring. For a Normal vertex take the next vertex which can be either an intersection
+# or an edge end.
+getnextresultvertex(v::RingVertex{Normal}) = nextvertex(v)
+getnextresultvertex(v::RingVertex{Entering}) = getnextclippedvertex(v)
+getnextresultvertex(v::RingVertex{Exiting}) = getnextclippingvertex(v)
 
 function clip(poly::Polygon, ring::Ring, ::WeilerAthertonClipping)
   polyrings = rings(poly)
@@ -85,35 +115,39 @@ function clip(poly::Polygon, ring::Ring, ::WeilerAthertonClipping)
     # Three consecutive clipping vertices are used to properly identify intersection vertex type
     # for corner cases, so the clipping segment is constructed with the two following vertices
     # after the current one.
-    clippingsegment = Segment(clipping.left.point, clipping.left.left.point)
+    clippingedgestart = nextnormal(clipping)
+    clippingedgeend = nextnormal(clippingedgestart)
+    clippingsegment = Segment(clippingedgestart.point, clippingedgeend.point)
 
     for (k, startclipped) in enumerate(clippedrings)
       clipped = startclipped
       while true
         # Like for the clipping, the clipped also uses three consecutive vertices.
-        clippedsegment = Segment(clipped.left.point, clipped.left.left.point)
+        clippededgestart = nextnormal(clipped)
+        clippededgeend = nextnormal(clippededgestart)
+        clippedsegment = Segment(clippededgestart.point, clippededgeend.point)
 
         I = intersection(clippingsegment, clippedsegment)
         vertex = vertexfromintersection(I, clipping, clipped)
-        success = insertintersections!(vertex, clipping, clipped, entering)
+        success = insertintersections!(vertex, clippingedgestart, clippededgestart, entering)
         intersected[k] = intersected[k] || success
 
-        clipped = clipped.left
-        clipped.left == startclipped.left && break
+        clipped = nextnormal(clipped)
+        nextnormal(clipped) == nextnormal(startclipped) && break
       end
     end
 
-    clipping = clipping.left
-    clipping.left == startclipping.left && break
+    clipping = nextnormal(clipping)
+    nextnormal(clipping) == nextnormal(startclipping) && break
   end
 
-  # Handle the case when no interections have been found.
   if !any(intersected)
+    # Handle the case when no interections have been found.
     if orientation(ring) == CW
-      # For an inner clipping ring take all outside rings.
+      # For an inner clipping ring take all clipped rings outside of it.
       return PolyArea(collectoutsiderings(ring, polyrings)...)
     else
-      # For an outer clipping ring add it to act as the outer ring.
+      # For an outer clipping ring add it to act as the outer ring of the clipped ring.
       collectedrings = all(v ∈ PolyArea(polyrings[1]) for v in vertices(ring)) ? [ring] : []
     end
   else
@@ -177,34 +211,37 @@ function perhapsaddnonintersected!(ring, polyrings, intersected)
   newpolyrings
 end
 
-# Inserts the intersection in the ring.
-function insertintersection!(head::RingVertex, intersection::RingVertex, side::Symbol)
-  tail = head.left
-  vertex = head
-
-  new = measure(Segment(head.point, intersection.point))
+function insertintersection!(head::RingVertex{Normal}, intersection::RingVertex, getnext, setnext!)
+  tail = nextnormal(head)
+  current = head
+  newdistance = measure(Segment(head.point, intersection.point))
+  # Search for the place to insert the intersection by comparing its distance
+  # from the head with other inserted intersections.
   while true
-    os = isnormal(vertex) ? :right : side
-    current = measure(Segment(head.point, getfield(vertex, os).point))
-    if (new < current) || (getfield(vertex, os) == tail)
-      next = getfield(vertex, os)
-      if !(new ≈ measure(Segment(head.point, next.point)))
-        setfield!(intersection, side, next)
-        setfield!(vertex, os, intersection)
+    before = current
+    current = getnext(current)
+    currentdistance = measure(Segment(head.point, current.point))
+
+    if (newdistance < currentdistance) || (current == tail)
+      if !(newdistance ≈ currentdistance) # TODO Check if this is needed.
+        # Only add if it is not coincident with the following vertex, or it could be added twice,
+        # if it is the tail vertex. In such case, the intersection will be detected second time
+        # and added then.
+        setnext!(intersection, current)
+        setnext!(before, intersection)
       end
       break
     end
-    vertex = getfield(vertex, os)
   end
 end
 
 # Inserts the intersection into both the clipping and the clipped rings.
-function insertintersections!(vertex::Tuple, clipping, clipped, entering)
+function insertintersections!(vertex::Tuple, clipping::RingVertex{Normal}, clipped::RingVertex{Normal}, entering)
   (vtype, point) = vertex
   if !isnothing(vtype)
     intersection = RingVertex{vtype}(point)
-    insertintersection!(clipping.left, intersection, :left)
-    insertintersection!(clipped.left, intersection, :right)
+    insertintersection!(clipping, intersection, getnextclippingvertex, setnextclippingvertex!)
+    insertintersection!(clipped, intersection, getnextclippedvertex, setnextclippedvertex!)
 
     if vtype == Entering
       push!(entering, intersection)
@@ -235,7 +272,7 @@ function collectclipped(entering::Vector{RingVertex{Entering}})
 
       push!(ring, vertex)
       push!(visited, vertex)
-      vertex = nextvertex(vertex)
+      vertex = getnextresultvertex(vertex)
     end
 
     # Remove duplicates.
@@ -273,30 +310,30 @@ function vertexfromintersection(I, clipping, clipped)
 end
 
 function vertexfromcrossing(point, clipping, clipped)
-  cl = Line(clipping.left.point, clipping.left.left.point)
-  vertextype = sideof(clipped.left.left.point, cl) == LEFT ? Entering : Exiting
+  cl = Line(nextnormal(clipping).point, nextnormal(nextnormal(clipping)).point)
+  vertextype = sideof(nextnormal(nextnormal(clipped)).point, cl) == LEFT ? Entering : Exiting
   (vertextype, point)
 end
 
 function vertexfromedgetouching(point, clipping, clipped)
   vertextype = nothing
-  if point ≈ clipped.left.point
+  if point ≈ nextnormal(clipped).point
     # When intersection is at the shared vertex of two edges of the clipped ring,
-    # then split the interscting edge of the clipping ring at the intersection point.
+    # then split the intersecting edge of the clipping ring at the intersection point.
     vertextype = decidedirection(
-      Segment(clipped.point, clipped.left.point),
-      Segment(clipped.left.point, clipped.left.left.point),
-      Segment(clipping.left.point, point),
-      Segment(point, clipping.left.left.point)
+      Segment(clipped.point, nextnormal(clipped).point),
+      Segment(nextnormal(clipped).point, nextnormal(nextnormal(clipped)).point),
+      Segment(nextnormal(clipping).point, point),
+      Segment(point, nextnormal(nextnormal(clipping)).point)
     )
-  elseif point ≈ clipping.left.point
+  elseif point ≈ nextnormal(clipping).point
     # When intersection is at the shared vertex of two edges of the clipping ring,
-    # then split the interscting edge of the clipped ring at the intersection point.
+    # then split the intersecting edge of the clipped ring at the intersection point.
     vertextype = decidedirection(
-      Segment(clipped.left.point, point),
-      Segment(point, clipped.left.left.point),
-      Segment(clipping.point, clipping.left.point),
-      Segment(clipping.left.point, clipping.left.left.point)
+      Segment(nextnormal(clipped).point, point),
+      Segment(point, nextnormal(nextnormal(clipped)).point),
+      Segment(clipping.point, nextnormal(clipping).point),
+      Segment(nextnormal(clipping).point, nextnormal(nextnormal(clipping)).point)
     )
   end
   (vertextype, point)
@@ -305,14 +342,14 @@ end
 function vertexfromcornertouching(point, clipping, clipped)
   vertextype = nothing
   # When intersection is at the shared vertices of both the clipping and the clipped rings.
-  if (point ≈ clipped.left.point) && (point ≈ clipping.left.point)
+  if (point ≈ nextnormal(clipped).point) && (point ≈ nextnormal(clipping).point)
     # Only applies if the intersection coincides with the middles of the currently observed
     # vertices for both clipping and clipped rings.
     vertextype = decidedirection(
       Segment(clipped.point, point),
-      Segment(point, clipped.left.left.point),
+      Segment(point, nextnormal(nextnormal(clipped)).point),
       Segment(clipping.point, point),
-      Segment(point, clipping.left.left.point)
+      Segment(point, nextnormal(nextnormal(clipping)).point)
     )
   end
   (vertextype, point)
@@ -320,35 +357,35 @@ end
 
 function vertexfromoverlapping(segment, clipping, clipped)
   # For both ends of the intersecting segment, check if it coincides with the middle of
-  # the obeserved vertices for clipped and clipping rings. If it does, attempt adding a
+  # the observed vertices for clipped and clipping rings. If it does, attempt adding a
   # point.
 
   ret = Tuple[]
   for point in extrema(segment)
-    if point ≈ clipped.left.point
-      clippingprev = Segment(clipping.left.point, point)
+    if point ≈ nextnormal(clipped).point
+      clippingprev = Segment(nextnormal(clipping).point, point)
       if measure(clippingprev) ≈ 0.0u"m"
         clippingprev = Segment(clipping.point, point)
       end
       vertextype = decidedirection(
         Segment(clipped.point, point),
-        Segment(point, clipped.left.left.point),
+        Segment(point, nextnormal(nextnormal(clipped)).point),
         clippingprev,
-        Segment(point, clipping.left.left.point)
+        Segment(point, nextnormal(nextnormal(clipping)).point)
       )
       push!(ret, (vertextype, point))
     end
 
-    if point ≈ clipping.left.point
-      clippedprev = Segment(clipped.left.point, point)
+    if point ≈ nextnormal(clipping).point
+      clippedprev = Segment(nextnormal(clipped).point, point)
       if measure(clippedprev) ≈ 0.0u"m"
         clippedprev = Segment(clipped.point, point)
       end
       vertextype = decidedirection(
         clippedprev,
-        Segment(point, clipped.left.left.point),
+        Segment(point, nextnormal(nextnormal(clipped)).point),
         Segment(clipping.point, point),
-        Segment(point, clipping.left.left.point)
+        Segment(point, nextnormal(nextnormal(clipping)).point)
       )
       push!(ret, (vertextype, point))
     end
@@ -416,5 +453,5 @@ function gettraversalring(ring::Ring)
     start = new
   end
 
-  return start.left
+  return nextnormal(start)
 end

From 538e0b3c3d2f2e4b0e476b88d7fb30a89974bfb1 Mon Sep 17 00:00:00 2001
From: Nemanja Jovanovic <nemanja.jovanovic@proton.me>
Date: Sun, 8 Dec 2024 18:09:10 +0200
Subject: [PATCH 5/6] Refactor intersection inserting loop Improve clipping in
 clipped test Improve points deduplication and removal of repeated entering
 vertices

---
 src/clipping/weileratherton.jl | 119 ++++++++++++++++-----------------
 test/clipping.jl               |   6 +-
 2 files changed, 57 insertions(+), 68 deletions(-)

diff --git a/src/clipping/weileratherton.jl b/src/clipping/weileratherton.jl
index 1cf04a22a..4ae7a9d27 100644
--- a/src/clipping/weileratherton.jl
+++ b/src/clipping/weileratherton.jl
@@ -46,6 +46,9 @@ RingVertex{VT}(point::Point{M,C}) where {VT<:VertexType,M<:Manifold,C<:CRS} = Ri
 isnormal(::RingVertex{Normal}) = true
 isnormal(::RingVertex) = false
 
+isentering(::RingVertex{Entering}) = true
+isentering(::RingVertex) = false
+
 function appendvertices!(v₁::RingVertex{Normal}, v₂::RingVertex{Normal})
   v₂.left = v₁.left
   v₁.left = v₂
@@ -100,18 +103,20 @@ function clip(poly::Polygon, ring::Ring, ::WeilerAthertonClipping)
 
   # Convert the subject polygon rings and the clipping ring to the RingVertex data structure.
   clippedrings = [gettraversalring(r) for r in polyrings]
-  startclipping = gettraversalring(ring)
+  clipping = gettraversalring(ring)
 
   # For keeping track of intersected rings, as the non-intersected ones need additional
   # processing at the end.
   intersected = zeros(Bool, length(clippedrings))
 
+  # For marking of clipping ring vertices which touch rings of clipped polygon.
+  clippingtouches = zeros(Bool, nvertices(ring))
+
   # For collecting all entering vertices, as they are used as starting points for collection
   # of the output rings.
   entering = RingVertex{Entering}[]
 
-  clipping = startclipping
-  while true
+  for l in 1:nvertices(ring)
     # Three consecutive clipping vertices are used to properly identify intersection vertex type
     # for corner cases, so the clipping segment is constructed with the two following vertices
     # after the current one.
@@ -119,26 +124,30 @@ function clip(poly::Polygon, ring::Ring, ::WeilerAthertonClipping)
     clippingedgeend = nextnormal(clippingedgestart)
     clippingsegment = Segment(clippingedgestart.point, clippingedgeend.point)
 
-    for (k, startclipped) in enumerate(clippedrings)
-      clipped = startclipped
-      while true
+    for (k, clipped) in enumerate(clippedrings)
+      for _ in 1:nvertices(polyrings[k])
         # Like for the clipping, the clipped also uses three consecutive vertices.
         clippededgestart = nextnormal(clipped)
         clippededgeend = nextnormal(clippededgestart)
         clippedsegment = Segment(clippededgestart.point, clippededgeend.point)
 
         I = intersection(clippingsegment, clippedsegment)
-        vertex = vertexfromintersection(I, clipping, clipped)
-        success = insertintersections!(vertex, clippingedgestart, clippededgestart, entering)
-        intersected[k] = intersected[k] || success
+
+        if type(I) != NotIntersecting
+          vertices = vertexfromintersection(I, clipping, clipped)
+          if !isnothing(vertices)
+            insertintersections!(vertices, clippingedgestart, clippededgestart, entering)
+            intersected[k] = true
+          elseif type(I) in [EdgeTouching, CornerTouching] && get(I) ≈ clippingedgestart.point
+            clippingtouches[l] = true
+          end
+        end
 
         clipped = nextnormal(clipped)
-        nextnormal(clipped) == nextnormal(startclipped) && break
       end
     end
 
     clipping = nextnormal(clipping)
-    nextnormal(clipping) == nextnormal(startclipping) && break
   end
 
   if !any(intersected)
@@ -147,8 +156,11 @@ function clip(poly::Polygon, ring::Ring, ::WeilerAthertonClipping)
       # For an inner clipping ring take all clipped rings outside of it.
       return PolyArea(collectoutsiderings(ring, polyrings)...)
     else
+      # Align the clippingtouches with the ring vertices.
+      clippingtouches = [clippingtouches[end], clippingtouches[1:(end - 1)]...]
       # For an outer clipping ring add it to act as the outer ring of the clipped ring.
-      collectedrings = all(v ∈ PolyArea(polyrings[1]) for v in vertices(ring)) ? [ring] : []
+      collectedrings =
+        all(t || v ∈ PolyArea(polyrings[1]) for (t, v) in zip(clippingtouches, vertices(ring))) ? [ring] : []
     end
   else
     # Collect rings formed from the intersected rings.
@@ -223,36 +235,27 @@ function insertintersection!(head::RingVertex{Normal}, intersection::RingVertex,
     currentdistance = measure(Segment(head.point, current.point))
 
     if (newdistance < currentdistance) || (current == tail)
-      if !(newdistance ≈ currentdistance) # TODO Check if this is needed.
-        # Only add if it is not coincident with the following vertex, or it could be added twice,
-        # if it is the tail vertex. In such case, the intersection will be detected second time
-        # and added then.
-        setnext!(intersection, current)
-        setnext!(before, intersection)
-      end
+      setnext!(intersection, current)
+      setnext!(before, intersection)
       break
     end
   end
 end
 
+insertintersections!(_::Nothing, _, _, _) = nothing
+
 # Inserts the intersection into both the clipping and the clipped rings.
-function insertintersections!(vertex::Tuple, clipping::RingVertex{Normal}, clipped::RingVertex{Normal}, entering)
-  (vtype, point) = vertex
-  if !isnothing(vtype)
-    intersection = RingVertex{vtype}(point)
-    insertintersection!(clipping, intersection, getnextclippingvertex, setnextclippingvertex!)
-    insertintersection!(clipped, intersection, getnextclippedvertex, setnextclippedvertex!)
-
-    if vtype == Entering
-      push!(entering, intersection)
-    end
-    return true
+function insertintersections!(vertex::RingVertex, clipping::RingVertex{Normal}, clipped::RingVertex{Normal}, entering)
+  insertintersection!(clipping, vertex, getnextclippingvertex, setnextclippingvertex!)
+  insertintersection!(clipped, vertex, getnextclippedvertex, setnextclippedvertex!)
+
+  if isentering(vertex)
+    push!(entering, vertex)
   end
-  false
 end
 
 insertintersections!(vertices::Array, clipping, clipped, entering) =
-  any(insertintersections!.(vertices, Ref(clipping), Ref(clipped), Ref(entering)))
+  insertintersections!.(vertices, Ref(clipping), Ref(clipped), Ref(entering))
 
 # Takes a list of entering vertices and returns all rings that contain those vertices.
 function collectclipped(entering::Vector{RingVertex{Entering}})
@@ -270,28 +273,14 @@ function collectclipped(entering::Vector{RingVertex{Entering}})
         break
       end
 
-      push!(ring, vertex)
-      push!(visited, vertex)
-      vertex = getnextresultvertex(vertex)
-    end
-
-    # Remove duplicates.
-    newring = RingVertex[ring[1]]
-    for i in 2:length(ring)
-      if !(ring[i].point ≈ newring[end].point)
-        push!(newring, ring[i])
+      nextvertex = getnextresultvertex(vertex)
+      # Skip adding sequential duplicates, and sequential entering vertices
+      # which can happen with overlapping edge intersections.
+      if !(isentering(vertex) && isentering(nextvertex)) && !(vertex.point ≈ nextvertex.point)
+        push!(ring, vertex)
       end
-    end
-    ring = newring
-
-    # Polygon might start several vertices after the first collected.
-    # This generally happens when there are overlapping edges that lead
-    # to several entering vertices without exiting in between. Then, the
-    # actual polygon is found by discarding the extra vertices before the
-    # proper loop.
-    k = findfirst(x -> ring[end].point == x.point, ring[1:(end - 1)])
-    if !isnothing(k)
-      ring = ring[(k + 1):end]
+      push!(visited, vertex)
+      vertex = nextvertex
     end
 
     if length(ring) > 2
@@ -306,13 +295,13 @@ function vertexfromintersection(I, clipping, clipped)
   type(I) == CornerTouching && return vertexfromcornertouching(get(I), clipping, clipped)
   type(I) == EdgeTouching && return vertexfromedgetouching(get(I), clipping, clipped)
   type(I) == Overlapping && return vertexfromoverlapping(get(I), clipping, clipped)
-  (nothing, nothing)
+  @assert false # No other intersection types expected.
 end
 
 function vertexfromcrossing(point, clipping, clipped)
   cl = Line(nextnormal(clipping).point, nextnormal(nextnormal(clipping)).point)
   vertextype = sideof(nextnormal(nextnormal(clipped)).point, cl) == LEFT ? Entering : Exiting
-  (vertextype, point)
+  RingVertex{vertextype}(point)
 end
 
 function vertexfromedgetouching(point, clipping, clipped)
@@ -336,7 +325,7 @@ function vertexfromedgetouching(point, clipping, clipped)
       Segment(nextnormal(clipping).point, nextnormal(nextnormal(clipping)).point)
     )
   end
-  (vertextype, point)
+  isnothing(vertextype) ? nothing : RingVertex{vertextype}(point)
 end
 
 function vertexfromcornertouching(point, clipping, clipped)
@@ -352,7 +341,7 @@ function vertexfromcornertouching(point, clipping, clipped)
       Segment(point, nextnormal(nextnormal(clipping)).point)
     )
   end
-  (vertextype, point)
+  isnothing(vertextype) ? nothing : RingVertex{vertextype}(point)
 end
 
 function vertexfromoverlapping(segment, clipping, clipped)
@@ -360,7 +349,7 @@ function vertexfromoverlapping(segment, clipping, clipped)
   # the observed vertices for clipped and clipping rings. If it does, attempt adding a
   # point.
 
-  ret = Tuple[]
+  ret = RingVertex[]
   for point in extrema(segment)
     if point ≈ nextnormal(clipped).point
       clippingprev = Segment(nextnormal(clipping).point, point)
@@ -373,7 +362,9 @@ function vertexfromoverlapping(segment, clipping, clipped)
         clippingprev,
         Segment(point, nextnormal(nextnormal(clipping)).point)
       )
-      push!(ret, (vertextype, point))
+      if !isnothing(vertextype)
+        push!(ret, RingVertex{vertextype}(point))
+      end
     end
 
     if point ≈ nextnormal(clipping).point
@@ -387,10 +378,12 @@ function vertexfromoverlapping(segment, clipping, clipped)
         Segment(clipping.point, point),
         Segment(point, nextnormal(nextnormal(clipping)).point)
       )
-      push!(ret, (vertextype, point))
+      if !isnothing(vertextype)
+        push!(ret, RingVertex{vertextype}(point))
+      end
     end
   end
-  ret
+  (length(ret) == 0) ? nothing : ret
 end
 
 # Used to figure out the type of the vertex to add when intersection is other than the crossing
@@ -419,13 +412,13 @@ function decidedirection(clipped₁, clipped₂, clipping₁, clipping₂)
   γ = mod(∠(a, c), twoπ)
   δ = mod(∠(a, d), twoπ)
 
-  if isapprox(γ, zero(γ), atol=tol) || isapprox(γ, twoπ, atol=tol)
+  if isapproxzero(γ, atol=tol) || isapprox(γ, twoπ, atol=tol)
     if δ < β < twoπ
       return Entering
     else
       return Exiting
     end
-  elseif isapprox(δ, zero(δ), atol=tol) || isapprox(δ, twoπ, atol=tol)
+  elseif isapproxzero(δ, atol=tol) || isapprox(δ, twoπ, atol=tol)
     if γ < β < twoπ
       return Exiting
     else
diff --git a/test/clipping.jl b/test/clipping.jl
index 7206afb75..34e6fd211 100644
--- a/test/clipping.jl
+++ b/test/clipping.jl
@@ -143,7 +143,7 @@ end
   clipped = clip(poly, other, WeilerAthertonClipping())
   @test !issimple(clipped)
   crings = rings(clipped)
-  @test all(vertices(crings[1]) .≈ [cart(4, 0), cart(4, 4), cart(0, 4), cart(0, 0)])
+  @test crings[1] ≗ Ring([cart(4, 0), cart(4, 4), cart(0, 4), cart(0, 0)])
   @test crings[2] ≗ inner
 
   # PolyArea with one inner ring inside `other` and another inner ring outside `other`
@@ -170,9 +170,6 @@ end
     ]
   )
 
-  # Tolerances are not properly retrieved for Float32 types, so need to pass them explicitly.
-  atol_el = coords(cart(0.0)).x
-
   # Overlapping clipping polygon with the edges of the hole of the subject polygon.
   triangle = Triangle((0.0, 0.0), (1.0, 0.0), (0.5, 1.0))
   rectangle = Quadrangle((0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0))
@@ -185,7 +182,6 @@ end
     isapprox.(
       vertices(crings[1]),
       [cart(0.8, 0.3), cart(0.3, 0.3), cart(0.15, 0.0), cart(0.95, 0.0)],
-      atol=atol(atol_el)
     )
   )
 end

From fa262641f6a5536ed535af459062bbc424fc9341 Mon Sep 17 00:00:00 2001
From: Nemanja Jovanovic <nemanja.jovanovic@proton.me>
Date: Fri, 17 Jan 2025 12:11:35 +0200
Subject: [PATCH 6/6] Simplify RingVertex point type

---
 src/clipping/weileratherton.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/clipping/weileratherton.jl b/src/clipping/weileratherton.jl
index 4ae7a9d27..9aa1ba517 100644
--- a/src/clipping/weileratherton.jl
+++ b/src/clipping/weileratherton.jl
@@ -28,12 +28,12 @@ abstract type Exiting <: VertexType end
 # Data structure for clipping the polygons. Fields left and right are used depending on the
 # VertexType, as designated with the helper functions. The data structure forms a directed
 # graph with each element always pointing to two elements.
-mutable struct RingVertex{VT<:VertexType,M<:Manifold,C<:CRS}
-  point::Point{M,C}
-  left::RingVertex{<:VertexType,M,C}
-  right::RingVertex{<:VertexType,M,C}
+mutable struct RingVertex{VT<:VertexType,P<:Point}
+  point::P
+  left::RingVertex{<:VertexType,P}
+  right::RingVertex{<:VertexType,P}
 
-  function RingVertex{VT,M,C}(point) where {VT<:VertexType,M<:Manifold,C<:CRS}
+  function RingVertex{VT,P}(point) where {VT<:VertexType,P<:Point}
     v = new(point)
     v.left = v
     v.right = v
@@ -41,7 +41,7 @@ mutable struct RingVertex{VT<:VertexType,M<:Manifold,C<:CRS}
   end
 end
 
-RingVertex{VT}(point::Point{M,C}) where {VT<:VertexType,M<:Manifold,C<:CRS} = RingVertex{VT,M,C}(point)
+RingVertex{VT}(point::P) where {VT<:VertexType,P<:Point} = RingVertex{VT,P}(point)
 
 isnormal(::RingVertex{Normal}) = true
 isnormal(::RingVertex) = false