Skip to content

Commit b09afad

Browse files
asinghvi17claude
asinghvi17
and
claude
authored
Experimental: Add Sutherland-Hodgman convex-convex polygon intersection (#375)
Co-authored-by: Claude Opus 4.5 <[email protected]>
1 parent d7d18ee commit b09afad

File tree

4 files changed

+339
-0
lines changed

4 files changed

+339
-0
lines changed

‎src/GeometryOps.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -74,6 +74,7 @@ include("methods/clipping/cut.jl")
7474
include("methods/clipping/intersection.jl")
7575
include("methods/clipping/difference.jl")
7676
include("methods/clipping/union.jl")
77+
include("methods/clipping/sutherland_hodgman.jl")
7778
include("methods/geom_relations/contains.jl")
7879
include("methods/geom_relations/coveredby.jl")
7980
include("methods/geom_relations/covers.jl")

‎src/methods/clipping/sutherland_hodgman.jl‎

Lines changed: 154 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,154 @@
1+
# # Sutherland-Hodgman Convex-Convex Clipping
2+
export ConvexConvexSutherlandHodgman
3+
4+
"""
5+
ConvexConvexSutherlandHodgman{M <: Manifold} <: GeometryOpsCore.Algorithm{M}
6+
7+
Sutherland-Hodgman polygon clipping algorithm optimized for convex-convex intersection.
8+
9+
Both input polygons MUST be convex. If either polygon is non-convex, results are undefined.
10+
11+
This is simpler and faster than Foster-Hormann for small convex polygons, with O(n*m)
12+
complexity where n and m are vertex counts.
13+
14+
# Example
15+
16+
```julia
17+
import GeometryOps as GO, GeoInterface as GI
18+
19+
square1 = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]])
20+
square2 = GI.Polygon([[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]])
21+
22+
result = GO.intersection(GO.ConvexConvexSutherlandHodgman(), square1, square2)
23+
```
24+
"""
25+
struct ConvexConvexSutherlandHodgman{M <: Manifold} <: GeometryOpsCore.Algorithm{M}
26+
manifold::M
27+
end
28+
29+
# Default constructor uses Planar
30+
ConvexConvexSutherlandHodgman() = ConvexConvexSutherlandHodgman(Planar())
31+
32+
# Main entry point - algorithm dispatch
33+
function intersection(
34+
alg::ConvexConvexSutherlandHodgman,
35+
geom_a,
36+
geom_b,
37+
::Type{T}=Float64;
38+
kwargs...
39+
) where {T<:AbstractFloat}
40+
return _intersection_sutherland_hodgman(
41+
alg, T,
42+
GI.trait(geom_a), geom_a,
43+
GI.trait(geom_b), geom_b
44+
)
45+
end
46+
47+
# Polygon-Polygon intersection using Sutherland-Hodgman
48+
function _intersection_sutherland_hodgman(
49+
alg::ConvexConvexSutherlandHodgman{Planar},
50+
::Type{T},
51+
::GI.PolygonTrait, poly_a,
52+
::GI.PolygonTrait, poly_b
53+
) where {T}
54+
# Get exterior rings (convex polygons have no holes)
55+
ring_a = GI.getexterior(poly_a)
56+
ring_b = GI.getexterior(poly_b)
57+
58+
# Start with vertices of poly_a as the output list (excluding closing point)
59+
output = Tuple{T,T}[]
60+
for point in GI.getpoint(ring_a)
61+
pt = _tuple_point(point, T)
62+
# Skip the closing point (same as first)
63+
if !isempty(output) && pt == output[1]
64+
continue
65+
end
66+
push!(output, pt)
67+
end
68+
69+
# Clip against each edge of poly_b
70+
for (edge_start, edge_end) in eachedge(ring_b, T)
71+
isempty(output) && break
72+
output = _sh_clip_to_edge(output, edge_start, edge_end, T)
73+
end
74+
75+
# Handle empty result (no intersection) - return degenerate polygon with zero area
76+
if isempty(output)
77+
zero_pt = (zero(T), zero(T))
78+
return GI.Polygon([[zero_pt, zero_pt, zero_pt]])
79+
end
80+
81+
# Close the ring
82+
push!(output, output[1])
83+
84+
# Return polygon
85+
return GI.Polygon([output])
86+
end
87+
88+
# Clip polygon against a single edge using Sutherland-Hodgman rules
89+
function _sh_clip_to_edge(polygon_points::Vector{Tuple{T,T}}, edge_start, edge_end, ::Type{T}) where T
90+
output = Tuple{T,T}[]
91+
n = length(polygon_points)
92+
n == 0 && return output
93+
94+
for i in 1:n
95+
current = polygon_points[i]
96+
next_pt = polygon_points[mod1(i + 1, n)]
97+
98+
# Determine if points are inside (left of or on the edge)
99+
# orient > 0 means left (inside for CCW polygon), == 0 means on edge, < 0 means right (outside)
100+
current_inside = Predicates.orient(edge_start, edge_end, current; exact=False()) >= 0
101+
next_inside = Predicates.orient(edge_start, edge_end, next_pt; exact=False()) >= 0
102+
103+
if current_inside
104+
push!(output, current)
105+
if !next_inside
106+
# Exiting: add intersection point
107+
intr_pt = _sh_line_intersection(current, next_pt, edge_start, edge_end, T)
108+
push!(output, intr_pt)
109+
end
110+
elseif next_inside
111+
# Entering: add intersection point
112+
intr_pt = _sh_line_intersection(current, next_pt, edge_start, edge_end, T)
113+
push!(output, intr_pt)
114+
end
115+
# Both outside: add nothing
116+
end
117+
118+
return output
119+
end
120+
121+
# Compute intersection point of line segment (p1, p2) with line through (p3, p4)
122+
function _sh_line_intersection(p1::Tuple{T,T}, p2::Tuple{T,T}, p3::Tuple{T,T}, p4::Tuple{T,T}, ::Type{T}) where T
123+
x1, y1 = p1
124+
x2, y2 = p2
125+
x3, y3 = p3
126+
x4, y4 = p4
127+
128+
denom = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4)
129+
130+
# Lines are parallel - shouldn't happen in valid Sutherland-Hodgman usage
131+
if abs(denom) < eps(T)
132+
return p1 # Fallback
133+
end
134+
135+
t = ((x1 - x3) * (y3 - y4) - (y1 - y3) * (x3 - x4)) / denom
136+
137+
x = x1 + t * (x2 - x1)
138+
y = y1 + t * (y2 - y1)
139+
140+
return (T(x), T(y))
141+
end
142+
143+
# Fallback for unsupported geometry combinations
144+
function _intersection_sutherland_hodgman(
145+
alg::ConvexConvexSutherlandHodgman,
146+
::Type{T},
147+
trait_a, geom_a,
148+
trait_b, geom_b
149+
) where {T}
150+
throw(ArgumentError(
151+
"ConvexConvexSutherlandHodgman only supports Polygon-Polygon intersection, " *
152+
"got $(typeof(trait_a)) and $(typeof(trait_b))"
153+
))
154+
end

‎test/methods/clipping/sutherland_hodgman.jl‎

Lines changed: 183 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,183 @@
1+
using Test
2+
import GeometryOps as GO
3+
import GeoInterface as GI
4+
5+
@testset "ConvexConvexSutherlandHodgman" begin
6+
@testset "Basic intersection" begin
7+
# Two overlapping squares - intersection is 1x1 square
8+
square1 = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]])
9+
square2 = GI.Polygon([[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]])
10+
11+
result = GO.intersection(GO.ConvexConvexSutherlandHodgman(), square1, square2)
12+
@test GI.trait(result) isa GI.PolygonTrait
13+
@test GO.area(result) ≈ 1.0 atol=1e-10
14+
end
15+
16+
@testset "No intersection" begin
17+
# Disjoint squares
18+
square1 = GI.Polygon([[(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)]])
19+
square2 = GI.Polygon([[(5.0, 5.0), (6.0, 5.0), (6.0, 6.0), (5.0, 6.0), (5.0, 5.0)]])
20+
21+
result = GO.intersection(GO.ConvexConvexSutherlandHodgman(), square1, square2)
22+
@test GI.trait(result) isa GI.PolygonTrait
23+
@test GO.area(result) ≈ 0.0 atol=1e-10
24+
end
25+
26+
@testset "One contains other" begin
27+
# Large square contains small square
28+
large = GI.Polygon([[(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)]])
29+
small = GI.Polygon([[(2.0, 2.0), (4.0, 2.0), (4.0, 4.0), (2.0, 4.0), (2.0, 2.0)]])
30+
31+
result = GO.intersection(GO.ConvexConvexSutherlandHodgman(), large, small)
32+
@test GI.trait(result) isa GI.PolygonTrait
33+
@test GO.area(result) ≈ 4.0 atol=1e-10
34+
35+
# Reverse order should give same result
36+
result2 = GO.intersection(GO.ConvexConvexSutherlandHodgman(), small, large)
37+
@test GO.area(result2) ≈ 4.0 atol=1e-10
38+
end
39+
40+
@testset "Triangles" begin
41+
# Two overlapping triangles (both CCW winding)
42+
tri1 = GI.Polygon([[(0.0, 0.0), (4.0, 0.0), (2.0, 4.0), (0.0, 0.0)]])
43+
tri2 = GI.Polygon([[(0.0, 2.0), (2.0, -2.0), (4.0, 2.0), (0.0, 2.0)]])
44+
45+
result = GO.intersection(GO.ConvexConvexSutherlandHodgman(), tri1, tri2)
46+
@test GI.trait(result) isa GI.PolygonTrait
47+
@test GO.area(result) > 0
48+
end
49+
50+
@testset "Identical polygons" begin
51+
# Same polygon should return itself
52+
square = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]])
53+
54+
result = GO.intersection(GO.ConvexConvexSutherlandHodgman(), square, square)
55+
@test GI.trait(result) isa GI.PolygonTrait
56+
@test GO.area(result) ≈ 4.0 atol=1e-10
57+
end
58+
59+
@testset "Shared edge" begin
60+
# Two squares sharing an edge
61+
square1 = GI.Polygon([[(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)]])
62+
square2 = GI.Polygon([[(1.0, 0.0), (2.0, 0.0), (2.0, 1.0), (1.0, 1.0), (1.0, 0.0)]])
63+
64+
result = GO.intersection(GO.ConvexConvexSutherlandHodgman(), square1, square2)
65+
@test GI.trait(result) isa GI.PolygonTrait
66+
# Shared edge only - area should be 0 or near 0
67+
@test GO.area(result) ≈ 0.0 atol=1e-10
68+
end
69+
70+
@testset "Unsupported geometry types" begin
71+
square = GI.Polygon([[(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)]])
72+
point = GI.Point(0.5, 0.5)
73+
74+
@test_throws ArgumentError GO.intersection(GO.ConvexConvexSutherlandHodgman(), square, point)
75+
@test_throws ArgumentError GO.intersection(GO.ConvexConvexSutherlandHodgman(), point, square)
76+
end
77+
78+
@testset "Pentagons" begin
79+
# Helper to create regular polygon (CCW winding)
80+
function regular_polygon(n, radius, center_x, center_y)
81+
coords = Tuple{Float64,Float64}[]
82+
for i in 0:n-1
83+
θ = 2π * i / n
84+
push!(coords, (center_x + radius * cos(θ), center_y + radius * sin(θ)))
85+
end
86+
push!(coords, coords[1]) # close the ring
87+
return GI.Polygon([coords])
88+
end
89+
90+
# Two overlapping pentagons
91+
pent1 = regular_polygon(5, 2.0, 0.0, 0.0)
92+
pent2 = regular_polygon(5, 2.0, 1.5, 0.0)
93+
94+
result = GO.intersection(GO.ConvexConvexSutherlandHodgman(), pent1, pent2)
95+
@test GI.trait(result) isa GI.PolygonTrait
96+
@test GO.area(result) > 0
97+
# Intersection should be smaller than either pentagon
98+
@test GO.area(result) < GO.area(pent1)
99+
end
100+
101+
@testset "Hexagons" begin
102+
function regular_polygon(n, radius, center_x, center_y)
103+
coords = Tuple{Float64,Float64}[]
104+
for i in 0:n-1
105+
θ = 2π * i / n
106+
push!(coords, (center_x + radius * cos(θ), center_y + radius * sin(θ)))
107+
end
108+
push!(coords, coords[1])
109+
return GI.Polygon([coords])
110+
end
111+
112+
# Two overlapping hexagons
113+
hex1 = regular_polygon(6, 2.0, 0.0, 0.0)
114+
hex2 = regular_polygon(6, 2.0, 2.0, 0.0)
115+
116+
result = GO.intersection(GO.ConvexConvexSutherlandHodgman(), hex1, hex2)
117+
@test GI.trait(result) isa GI.PolygonTrait
118+
@test GO.area(result) > 0
119+
@test GO.area(result) < GO.area(hex1)
120+
121+
# Hexagon containing smaller hexagon
122+
hex_large = regular_polygon(6, 3.0, 0.0, 0.0)
123+
hex_small = regular_polygon(6, 1.0, 0.0, 0.0)
124+
125+
result2 = GO.intersection(GO.ConvexConvexSutherlandHodgman(), hex_large, hex_small)
126+
@test GO.area(result2) ≈ GO.area(hex_small) atol=1e-10
127+
end
128+
129+
@testset "Octagons" begin
130+
function regular_polygon(n, radius, center_x, center_y)
131+
coords = Tuple{Float64,Float64}[]
132+
for i in 0:n-1
133+
θ = 2π * i / n
134+
push!(coords, (center_x + radius * cos(θ), center_y + radius * sin(θ)))
135+
end
136+
push!(coords, coords[1])
137+
return GI.Polygon([coords])
138+
end
139+
140+
# Two overlapping octagons
141+
oct1 = regular_polygon(8, 2.0, 0.0, 0.0)
142+
oct2 = regular_polygon(8, 2.0, 1.0, 1.0)
143+
144+
result = GO.intersection(GO.ConvexConvexSutherlandHodgman(), oct1, oct2)
145+
@test GI.trait(result) isa GI.PolygonTrait
146+
@test GO.area(result) > 0
147+
@test GO.area(result) < GO.area(oct1)
148+
149+
# Identical octagons should return same area
150+
result2 = GO.intersection(GO.ConvexConvexSutherlandHodgman(), oct1, oct1)
151+
@test GO.area(result2) ≈ GO.area(oct1) atol=1e-10
152+
end
153+
154+
@testset "Bordering rectangles with offset" begin
155+
# Two rectangles sharing an edge but offset vertically
156+
# rect1: [0,2] x [0,2], rect2: [2,4] x [0.5,2.5]
157+
# They share the edge at x=2 but are offset by 0.5 in y
158+
# Intersection should be zero area (just a line segment)
159+
rect1 = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]])
160+
rect2 = GI.Polygon([[(2.0, 0.5), (4.0, 0.5), (4.0, 2.5), (2.0, 2.5), (2.0, 0.5)]])
161+
162+
result = GO.intersection(GO.ConvexConvexSutherlandHodgman(), rect1, rect2)
163+
@test GI.trait(result) isa GI.PolygonTrait
164+
@test GO.area(result) ≈ 0.0 atol=1e-10
165+
166+
# Two rectangles sharing an edge but offset horizontally
167+
# rect3: [0,2] x [0,2], rect4: [0.5,2.5] x [2,4]
168+
rect3 = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]])
169+
rect4 = GI.Polygon([[(0.5, 2.0), (2.5, 2.0), (2.5, 4.0), (0.5, 4.0), (0.5, 2.0)]])
170+
171+
result2 = GO.intersection(GO.ConvexConvexSutherlandHodgman(), rect3, rect4)
172+
@test GI.trait(result2) isa GI.PolygonTrait
173+
@test GO.area(result2) ≈ 0.0 atol=1e-10
174+
175+
# Rectangles that just touch at a corner
176+
rect5 = GI.Polygon([[(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)]])
177+
rect6 = GI.Polygon([[(1.0, 1.0), (2.0, 1.0), (2.0, 2.0), (1.0, 2.0), (1.0, 1.0)]])
178+
179+
result3 = GO.intersection(GO.ConvexConvexSutherlandHodgman(), rect5, rect6)
180+
@test GI.trait(result3) isa GI.PolygonTrait
181+
@test GO.area(result3) ≈ 0.0 atol=1e-10
182+
end
183+
end

‎test/runtests.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,7 @@ end
3838
@safetestset "Cut" begin include("methods/clipping/cut.jl") end
3939
@safetestset "Intersection Point" begin include("methods/clipping/intersection_points.jl") end
4040
@safetestset "Polygon Clipping" begin include("methods/clipping/polygon_clipping.jl") end
41+
@safetestset "Sutherland-Hodgman" begin include("methods/clipping/sutherland_hodgman.jl") end
4142
# Transformations
4243
@safetestset "Embed Extent" begin include("transformations/extent.jl") end
4344
@safetestset "Reproject" begin include("transformations/reproject.jl") end

0 commit comments

Comments
 (0)