@@ -359,296 +359,3 @@ function _intersection_sutherland_hodgman(
359359 " got $(typeof (trait_a)) and $(typeof (trait_b)) "
360360 ))
361361end
362-
363- # =============================================================================
364- # Spherical Implementation
365- # =============================================================================
366-
367- using .. UnitSpherical: UnitSphericalPoint, UnitSphereFromGeographic, GeographicFromUnitSphere
368- using LinearAlgebra: cross, dot, normalize, norm
369-
370- """
371- spherical_orient(a, b, c)
372-
373- Compute the orientation of point `c` with respect to the great circle through `a` and `b`.
374-
375- Returns:
376- - Positive if `c` is to the left of the directed great circle arc from `a` to `b`
377- - Negative if `c` is to the right
378- - Zero if `c` is on the great circle
379-
380- The orientation is computed as the triple scalar product: (a × b) · c
381- """
382- function spherical_orient (a:: UnitSphericalPoint , b:: UnitSphericalPoint , c:: UnitSphericalPoint )
383- # Triple scalar product: (a × b) · c
384- return dot (cross (a, b), c)
385- end
386-
387- """
388- _point_in_convex_spherical_polygon(point, polygon_points)
389-
390- Check if a point is STRICTLY inside a convex spherical polygon.
391-
392- Returns `true` only if the point is strictly inside (all orientations positive).
393- Returns `false` if the point is on the boundary (any orientation is zero) or outside.
394-
395- This is critical for the fallback containment check in Sutherland-Hodgman:
396- points exactly ON the boundary should NOT be considered "inside", as this would
397- cause adjacent (non-overlapping) polygons to incorrectly return the clip polygon.
398- """
399- function _point_in_convex_spherical_polygon (
400- point:: UnitSphericalPoint ,
401- polygon_points:: Vector{<:UnitSphericalPoint}
402- )
403- n = length (polygon_points)
404- n < 3 && return false
405-
406- for i in 1 : n
407- edge_start = polygon_points[i]
408- edge_end = polygon_points[mod1 (i + 1 , n)]
409- orient = spherical_orient (edge_start, edge_end, point)
410- if orient <= 0 # strictly outside OR on the boundary
411- return false
412- end
413- end
414- return true # strictly inside (all orient > 0)
415- end
416-
417- # Helper to convert a GeoInterface point to UnitSphericalPoint
418- function _to_unit_spherical (point)
419- return UnitSphericalPoint (GI. PointTrait (), point)
420- end
421-
422- # Spherical Polygon-Polygon intersection using Sutherland-Hodgman
423- function _intersection_sutherland_hodgman (
424- alg:: ConvexConvexSutherlandHodgman{<:Spherical} ,
425- :: Type{T} ,
426- :: GI.PolygonTrait , poly_a,
427- :: GI.PolygonTrait , poly_b
428- ) where {T}
429- # Get exterior rings (convex polygons have no holes)
430- ring_a = GI. getexterior (poly_a)
431- ring_b = GI. getexterior (poly_b)
432-
433- # Convert poly_a vertices to UnitSphericalPoints (excluding closing point)
434- output = UnitSphericalPoint{T}[]
435- for point in GI. getpoint (ring_a)
436- pt = _to_unit_spherical (point)
437- pt_typed = UnitSphericalPoint {T} (pt. x, pt. y, pt. z)
438- # Skip the closing point (same as first)
439- if ! isempty (output) && isapprox (pt_typed, output[1 ]; atol= eps (T)* 10 )
440- continue
441- end
442- push! (output, pt_typed)
443- end
444-
445- # Store original subject polygon for containment check
446- original_subject = copy (output)
447-
448- # Convert clip polygon vertices (excluding closing point)
449- clip_points = UnitSphericalPoint{T}[]
450- for point in GI. getpoint (ring_b)
451- pt = _to_unit_spherical (point)
452- pt_typed = UnitSphericalPoint {T} (pt. x, pt. y, pt. z)
453- if ! isempty (clip_points) && isapprox (pt_typed, clip_points[1 ]; atol= eps (T)* 10 )
454- continue
455- end
456- push! (clip_points, pt_typed)
457- end
458-
459- # Clip against each edge of poly_b
460- n_clip = length (clip_points)
461- for i in 1 : n_clip
462- isempty (output) && break
463- edge_start = clip_points[i]
464- edge_end = clip_points[mod1 (i + 1 , n_clip)]
465- output = _sh_spherical_clip_to_edge (output, edge_start, edge_end, T)
466- end
467-
468- # Handle empty result - check if clip polygon is fully inside the original subject
469- if isempty (output)
470- # Use strict interior check: point ON boundary returns false
471- if ! isempty (clip_points) && _point_in_convex_spherical_polygon (clip_points[1 ], original_subject)
472- # Subject strictly contains clip - return clip polygon
473- result = copy (clip_points)
474- push! (result, result[1 ])
475- return _spherical_polygon_to_geo (result, T)
476- end
477- # No intersection - return degenerate polygon with zero area
478- zero_pt = (zero (T), zero (T))
479- return GI. Polygon ([[zero_pt, zero_pt, zero_pt]])
480- end
481-
482- # Check for degenerate result (collinear points = line segment, not polygon)
483- # This happens when adjacent polygons share an edge
484- if _is_degenerate_spherical_polygon (output, T)
485- zero_pt = (zero (T), zero (T))
486- return GI. Polygon ([[zero_pt, zero_pt, zero_pt]])
487- end
488-
489- # Close the ring and convert back to geographic coordinates
490- push! (output, output[1 ])
491- return _spherical_polygon_to_geo (output, T)
492- end
493-
494- """
495- _is_degenerate_spherical_polygon(points, T)
496-
497- Check if a set of spherical points forms a degenerate polygon (zero area).
498-
499- A polygon is degenerate if:
500- 1. It has fewer than 3 distinct points
501- 2. All points are collinear (lie on a single great circle)
502-
503- This is critical for detecting edge-sharing cases where Sutherland-Hodgman
504- produces a "polygon" that is actually just a line segment.
505- """
506- function _is_degenerate_spherical_polygon (
507- points:: Vector{UnitSphericalPoint{T}} ,
508- :: Type{T}
509- ) where T
510- n = length (points)
511- n < 3 && return true
512-
513- # Remove duplicate/very close points
514- unique_points = UnitSphericalPoint{T}[]
515- for p in points
516- is_dup = false
517- for existing in unique_points
518- if isapprox (p, existing; atol= eps (T)* 100 )
519- is_dup = true
520- break
521- end
522- end
523- if ! is_dup
524- push! (unique_points, p)
525- end
526- end
527-
528- length (unique_points) < 3 && return true
529-
530- # Check if all points are collinear on the sphere
531- # Points are collinear if they all lie on the same great circle
532- # This is true if for any three points A, B, C: (A × B) · C ≈ 0
533- p1 = unique_points[1 ]
534- p2 = unique_points[2 ]
535-
536- # Normal to the great circle through p1 and p2
537- normal = cross (p1, p2)
538- normal_len_sq = dot (normal, normal)
539-
540- # If p1 and p2 are nearly identical or antipodal, can't determine great circle
541- if normal_len_sq < eps (T)
542- return true
543- end
544-
545- # Check if all other points lie on the same great circle
546- for i in 3 : length (unique_points)
547- p = unique_points[i]
548- # Distance from point to great circle plane
549- dist = abs (dot (normal, p)) / sqrt (normal_len_sq)
550- if dist > sqrt (eps (T)) # Point is not on the great circle
551- return false # Non-degenerate polygon
552- end
553- end
554-
555- return true # All points are collinear
556- end
557-
558- # Convert UnitSphericalPoints back to geographic polygon
559- function _spherical_polygon_to_geo (points:: Vector{UnitSphericalPoint{T}} , :: Type{T} ) where T
560- geo_transform = GeographicFromUnitSphere ()
561- geo_points = [geo_transform (p) for p in points]
562- return GI. Polygon ([geo_points])
563- end
564-
565- # Clip polygon against a single spherical edge using Sutherland-Hodgman rules
566- function _sh_spherical_clip_to_edge (
567- polygon_points:: Vector{UnitSphericalPoint{T}} ,
568- edge_start:: UnitSphericalPoint{T} ,
569- edge_end:: UnitSphericalPoint{T} ,
570- :: Type{T}
571- ) where T
572- output = UnitSphericalPoint{T}[]
573- n = length (polygon_points)
574- n == 0 && return output
575-
576- for i in 1 : n
577- current = polygon_points[i]
578- next_pt = polygon_points[mod1 (i + 1 , n)]
579-
580- # Determine if points are inside (left of or on the edge)
581- # orient > 0 means left (inside for CCW polygon), == 0 means on edge, < 0 means right (outside)
582- current_orient = spherical_orient (edge_start, edge_end, current)
583- next_orient = spherical_orient (edge_start, edge_end, next_pt)
584-
585- current_inside = current_orient >= 0
586- next_inside = next_orient >= 0
587-
588- if current_inside
589- push! (output, current)
590- if ! next_inside
591- # Exiting: add intersection point
592- intr_pt = _sh_spherical_line_intersection (current, next_pt, edge_start, edge_end, T)
593- if ! isnothing (intr_pt)
594- push! (output, intr_pt)
595- end
596- end
597- elseif next_inside
598- # Entering: add intersection point
599- intr_pt = _sh_spherical_line_intersection (current, next_pt, edge_start, edge_end, T)
600- if ! isnothing (intr_pt)
601- push! (output, intr_pt)
602- end
603- end
604- # Both outside: add nothing
605- end
606-
607- return output
608- end
609-
610- # Compute intersection point of two great circle arcs
611- function _sh_spherical_line_intersection (
612- p1:: UnitSphericalPoint{T} ,
613- p2:: UnitSphericalPoint{T} ,
614- p3:: UnitSphericalPoint{T} ,
615- p4:: UnitSphericalPoint{T} ,
616- :: Type{T}
617- ) where T
618- # The intersection of two great circles is found by:
619- # 1. Compute normal vectors to each great circle plane
620- # 2. The intersection is the cross product of these normals (normalized)
621-
622- # Normal to great circle through p1, p2
623- n1 = cross (p1, p2)
624- # Normal to great circle through p3, p4
625- n2 = cross (p3, p4)
626-
627- # Cross product gives intersection direction
628- intersection_dir = cross (n1, n2)
629-
630- # Check if great circles are parallel (no intersection or same circle)
631- len_sq = dot (intersection_dir, intersection_dir)
632- if len_sq < eps (T)^ 2
633- # Great circles are parallel - return midpoint as fallback
634- return nothing
635- end
636-
637- # Normalize to get point on unit sphere
638- intersection_pt = normalize (intersection_dir)
639-
640- # The cross product gives us one of two antipodal points
641- # We need to pick the one that lies on the arcs (or closer to them)
642- # Check if this point is on the correct side of both arcs
643- # Use dot product to check if point is in the "positive" direction
644-
645- # Check which of the two antipodal points is closer to the arc midpoints
646- mid1 = normalize (p1 + p2)
647- mid2 = normalize (p3 + p4)
648-
649- if dot (intersection_pt, mid1) < 0 || dot (intersection_pt, mid2) < 0
650- intersection_pt = - intersection_pt
651- end
652-
653- return UnitSphericalPoint {T} (intersection_pt. x, intersection_pt. y, intersection_pt. z)
654- end
0 commit comments