Skip to content

Commit c899620

Browse files
committed
Extract shared Hao-Sun algorithm logic into helper function
Address PR #340 review comments by: - Creating `_hao_sun_edge_case` helper function to eliminate code duplication - Adding detailed comments explaining the 26 cases from Hao & Sun (2018) - Refactoring both standard and tree-accelerated implementations to use shared logic
1 parent 7677bf9 commit c899620

File tree

1 file changed

+62
-43
lines changed

1 file changed

+62
-43
lines changed

src/methods/geom_relations/geom_geom_processors.jl

Lines changed: 62 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -480,7 +480,7 @@ linestring or linearring trait, that is assumed to be closed, regardless of
480480
repeated last point.
481481
482482
Can provide values of in, on, and out keywords, which determines return values
483-
for each scenario.
483+
for each scenario.
484484
485485
Note that this uses the Algorithm by Hao and Sun (2018):
486486
https://doi.org/10.3390/sym10100477
@@ -493,6 +493,60 @@ of the curve if it didn't return 'on'.
493493
See paper for more information on cases denoted in comments.
494494
=#
495495

496+
#=
497+
Helper function that implements the core Hao-Sun algorithm logic for a single edge.
498+
This function processes one edge of a polygon to determine if a ray from point (x,y)
499+
intersects the edge, which is used for point-in-polygon testing.
500+
501+
Arguments:
502+
- x, y: coordinates of the test point
503+
- p_start, p_end: start and end points of the edge
504+
- exact: whether to use exact predicates
505+
506+
Returns a tuple (should_return_on, should_increment_k):
507+
- should_return_on: true if the point is on the edge (early return case)
508+
- should_increment_k: true if the ray crosses this edge (increment crossing counter)
509+
510+
The function implements the 26 cases from Hao & Sun (2018):
511+
- Cases 1-2: Point on horizontal edge
512+
- Cases 3-4, 9-10: Ray crosses edge (not at vertex)
513+
- Cases 7-8, 14-15, 17-18: Ray touches edge at vertex (various configurations)
514+
- Cases 11, 26: Ray doesn't intersect edge (both vertices on same side)
515+
- Cases 12-13, 16, 19-25: Ray intersects at edge vertex (need orientation check)
516+
=#
517+
@inline function _hao_sun_edge_case(x, y, p_start, p_end; exact)
518+
v1 = GI.y(p_start) - y
519+
v2 = GI.y(p_end) - y
520+
521+
# Skip if both vertices are on the same side of the horizontal ray (cases 11 or 26)
522+
((v1 < 0 && v2 < 0) || (v1 > 0 && v2 > 0)) && return (false, false)
523+
524+
u1, u2 = GI.x(p_start) - x, GI.x(p_end) - x
525+
f = Predicates.orient(p_start, p_end, (x, y); exact)
526+
527+
if v2 > 0 && v1 0 # Cases 3, 9, 16, 21, 13, or 24
528+
# Edge crosses ray upward
529+
f == 0 && return (true, false) # Point on edge (cases 16 or 21)
530+
f > 0 && return (false, true) # Ray crosses edge from right (cases 3 or 9)
531+
elseif v1 > 0 && v2 0 # Cases 4, 10, 19, 20, 12, or 25
532+
# Edge crosses ray downward
533+
f == 0 && return (true, false) # Point on edge (cases 19 or 20)
534+
f < 0 && return (false, true) # Ray crosses edge from left (cases 4 or 10)
535+
elseif v2 == 0 && v1 < 0 # Cases 7, 14, or 17
536+
# End vertex on ray, start vertex below
537+
f == 0 && return (true, false) # Point on edge (case 17)
538+
elseif v1 == 0 && v2 < 0 # Cases 8, 15, or 18
539+
# Start vertex on ray, end vertex below
540+
f == 0 && return (true, false) # Point on edge (case 18)
541+
elseif v1 == 0 && v2 == 0 # Cases 1, 2, 5, 6, 22, or 23
542+
# Both vertices on ray (horizontal edge)
543+
u2 0 && u1 0 && return (true, false) # Point between vertices (case 1)
544+
u1 0 && u2 0 && return (true, false) # Point between vertices (case 2)
545+
end
546+
547+
return (false, false) # No crossing, point not on edge
548+
end
549+
496550
function _point_filled_curve_orientation(
497551
::Planar, point, curve;
498552
in::T = point_in, on::T = point_on, out::T = point_out, exact,
@@ -504,26 +558,9 @@ function _point_filled_curve_orientation(
504558
p_start = GI.getpoint(curve, n)
505559
for (i, p_end) in enumerate(GI.getpoint(curve))
506560
i > n && break
507-
v1 = GI.y(p_start) - y
508-
v2 = GI.y(p_end) - y
509-
if !((v1 < 0 && v2 < 0) || (v1 > 0 && v2 > 0)) # if not cases 11 or 26
510-
u1, u2 = GI.x(p_start) - x, GI.x(p_end) - x
511-
f = Predicates.orient(p_start, p_end, (x, y); exact)
512-
if v2 > 0 && v1 0 # Case 3, 9, 16, 21, 13, or 24
513-
f == 0 && return on # Case 16 or 21
514-
f > 0 && (k += 1) # Case 3 or 9
515-
elseif v1 > 0 && v2 0 # Case 4, 10, 19, 20, 12, or 25
516-
f == 0 && return on # Case 19 or 20
517-
f < 0 && (k += 1) # Case 4 or 10
518-
elseif v2 == 0 && v1 < 0 # Case 7, 14, or 17
519-
f == 0 && return on # Case 17
520-
elseif v1 == 0 && v2 < 0 # Case 8, 15, or 18
521-
f == 0 && return on # Case 18
522-
elseif v1 == 0 && v2 == 0 # Case 1, 2, 5, 6, 22, or 23
523-
u2 0 && u1 0 && return on # Case 1
524-
u1 0 && u2 0 && return on # Case 2
525-
end
526-
end
561+
should_return_on, should_increment_k = _hao_sun_edge_case(x, y, p_start, p_end; exact)
562+
should_return_on && return on
563+
should_increment_k && (k += 1)
527564
p_start = p_end
528565
end
529566
return iseven(k) ? out : in
@@ -544,28 +581,10 @@ function _point_filled_curve_orientation(
544581
function per_edge_function(i)
545582
p_start = _tuple_point(GI.getpoint(curve, i))
546583
p_end = _tuple_point(GI.getpoint(curve, i + 1))
547-
v1 = GI.y(p_start) - y
548-
v2 = GI.y(p_end) - y
549-
if !((v1 < 0 && v2 < 0) || (v1 > 0 && v2 > 0)) # if not cases 11 or 26
550-
u1, u2 = GI.x(p_start) - x, GI.x(p_end) - x
551-
f = Predicates.orient(p_start, p_end, (x, y); exact)
552-
if v2 > 0 && v1 0 # Case 3, 9, 16, 21, 13, or 24
553-
f == 0 && return LSM.Action{T}(:full_return, on) # Case 16 or 21
554-
f > 0 && (k += 1) # Case 3 or 9
555-
elseif v1 > 0 && v2 0 # Case 4, 10, 19, 20, 12, or 25
556-
f == 0 && return LSM.Action{T}(:full_return, on) # Case 19 or 20
557-
f < 0 && (k += 1) # Case 4 or 10
558-
elseif v2 == 0 && v1 < 0 # Case 7, 14, or 17
559-
f == 0 && return LSM.Action{T}(:full_return, on) # Case 17
560-
elseif v1 == 0 && v2 < 0 # Case 8, 15, or 18
561-
f == 0 && return LSM.Action{T}(:full_return, on) # Case 18
562-
elseif v1 == 0 && v2 == 0 # Case 1, 2, 5, 6, 22, or 23
563-
u2 0 && u1 0 && return LSM.Action{T}(:full_return, on) # Case 1
564-
u1 0 && u2 0 && return LSM.Action{T}(:full_return, on) # Case 2
565-
end
566-
return LSM.Action(:continue, on)
567-
end
568-
p_start = p_end
584+
should_return_on, should_increment_k = _hao_sun_edge_case(x, y, p_start, p_end; exact)
585+
should_return_on && return LSM.Action{T}(:full_return, on)
586+
should_increment_k && (k += 1)
587+
return nothing
569588
end
570589

571590
result = SpatialTreeInterface.depth_first_search(per_edge_function,extent -> extent.Y[1] <= y <= extent.Y[2], tree)

0 commit comments

Comments
 (0)