Harden spherical routines for near-antipodal inputs (S2 parity) (#399)

Co-authored-by: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

Author: asinghvi17

src/utils/UnitSpherical/point.jl

@@ -114,8 +114,14 @@ GI.crs(::UnitSphericalPoint) = GFT.ProjString("+proj=cart +R=1 +type=crs") # TOD
 """
     spherical_distance(x::UnitSphericalPoint, y::UnitSphericalPoint)
 
-Compute the spherical distance between two points on the unit sphere.  
-Returns a `Number`, usually Float64 but that depends on the input type.
+Compute the great-circle distance (central angle, in radians) between two
+points on the unit sphere.  Returns a `Number` whose type follows from the
+inputs.
+
+Uses the `atan2(‖x × y‖, x · y)` form, which is numerically stable across
+the full `[0, π]` range — unlike `acos(x · y)`, which loses precision for
+nearly-identical points.  Adapted from Google's S2 geometry library
+(see [`Vector3::Angle`](https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/util/math/vector.h#L492)).
 
 # Extended help
 
@@ -135,7 +141,7 @@ julia> spherical_distance(UnitSphericalPoint(1, 0, 0), UnitSphericalPoint(1, 0,
 0.0
 ```
 """
-spherical_distance(x::UnitSphericalPoint, y::UnitSphericalPoint) = acos(clamp(x ⋅ y, -1.0, 1.0))
+spherical_distance(x::UnitSphericalPoint, y::UnitSphericalPoint) = atan(norm(cross(x, y)), x ⋅ y)
 
 # ## Random points
 Random.rand(rng::Random.AbstractRNG, ::Random.SamplerType{UnitSphericalPoint}) = rand(rng, UnitSphericalPoint{Float64})

src/utils/UnitSpherical/robustcrossproduct/RobustCrossProduct.jl

@@ -211,9 +211,13 @@ function exact_cross_product(a::AbstractVector{T}, b::AbstractVector{T}) where T
     big_b = BigFloat.(b; precision=512)
     result_xf = cross(big_a, big_b)
 
-    # Check if we got a non-zero result
-    # This is equivalent to s2's `s2pred::IsZero`.
-    if !all(<=(1e-300), abs.(result_xf))
+    # Check if we got a non-zero result. Matches S2's s2pred::IsZero at
+    # s2predicates_internal.h:77-79, which is a literal sign == 0 test on
+    # each component (i.e., all three components exactly zero). We must NOT
+    # use a magnitude threshold here; any non-zero BigFloat result — even
+    # deeply subnormal by Float64 standards — should go through the exact
+    # rescaling path in `normalizableFromExact`.
+    if !(iszero(result_xf[1]) && iszero(result_xf[2]) && iszero(result_xf[3]))
         return normalizableFromExact(result_xf)
     end
 
@@ -322,9 +326,11 @@ function symbolic_cross_product_sorted(a::AbstractVector{T}, b::AbstractVector{T
     @assert b[1] == 0 && b[2] == 0 "Expected both b[1] and b[2] to be zero"
 
     if a[1] != 0 || a[2] != 0
-        # Fix: This needs to match C++ code which returns (-a[1], a[0], 0) in 0-based indexing
-        # In Julia's 1-based indexing, this is (-a[2], a[1], 0)
-        return UnitSphericalPoint{T}(-a[2], a[1], 0)
+        # S2's SymbolicCrossProdSorted at s2edge_crossings.cc:251-253 returns
+        # `Vector3_d(a[1], -a[0], 0)` in C++ 0-based indexing.
+        # In Julia's 1-based indexing, that's (a[2], -a[1], 0). Bug fixed at
+        # src/utils/UnitSpherical/robustcrossproduct/RobustCrossProduct.jl:327.
+        return UnitSphericalPoint{T}(a[2], -a[1], 0)
     end
 
     # This is always non-zero in the S2 implementation

src/utils/UnitSpherical/robustcrossproduct/utils.jl

@@ -100,27 +100,45 @@ loss of precision due to floating-point underflow.
 This matches S2's NormalizableFromExact function.
 """
 function normalizableFromExact(xf::AbstractVector{BigFloat})
-    # First try a simple conversion
+    # Mirrors S2's NormalizableFromExact at s2edge_crossings.cc:318-336.
+    #
+    # First try a simple conversion to Float64; if that already has a
+    # component >= 2^-242 it needs no rescaling. Otherwise we must rescale
+    # *in BigFloat* (or equivalently compute the shift from BigFloat
+    # exponents) — we cannot convert to Float64 first, because components
+    # like 5e-324 would flush to the subnormal range or zero, destroying
+    # the axis information before the scaling multiply can fix it.
+    # (Bug fixed at src/utils/UnitSpherical/robustcrossproduct/utils.jl:102.)
     x = Float64.(xf)
-    
+
     if isNormalizable(x)
         return x
     end
-    
-    # Find the largest exponent
-    max_exp = -1000000  # Very small initial value
+
+    # Find the largest BigFloat exponent among nonzero components.
+    found = false
+    max_exp = 0
     for i in 1:3
         if !iszero(xf[i])
-            max_exp = max(max_exp, exponent(xf[i]))
+            e = exponent(xf[i])
+            if !found || e > max_exp
+                max_exp = e
+                found = true
+            end
         end
     end
-    
-    if max_exp < -1000000  # No non-zero components
-        return zero(xf)
+
+    if !found
+        return zero(x)  # The exact result is (0, 0, 0).
     end
-    
-    # Scale to get components in a good range
-    return Float64.(ldexp.(Float64.(xf), -max_exp))
+
+    # Scale in BigFloat so the largest component is in [0.5, 1), then
+    # convert. This matches S2's `ldexp(xf[i], -exp)` inside ExactFloat.
+    return UnitSphericalPoint(
+        Float64(ldexp(xf[1], -max_exp)),
+        Float64(ldexp(xf[2], -max_exp)),
+        Float64(ldexp(xf[3], -max_exp)),
+    )
 end
 
 

src/utils/UnitSpherical/slerp.jl

@@ -20,37 +20,53 @@ interpolates along that path.
 """
     slerp(a::UnitSphericalPoint, b::UnitSphericalPoint, i01::Number)
 
-Interpolate between `a` and `b`, at a proportion `i01` 
+Interpolate between `a` and `b`, at a proportion `i01`
 between 0 and 1 along the path from `a` to `b`.
 
+Uses the tangent-vector form `cos(r)·a + sin(r)·dir` — where `r = i01 ·
+spherical_distance(a, b)` and `dir = normalize(robust_cross_product(a, b) × a)`
+is the unit tangent at `a` pointing toward `b`. This avoids the `1/sin(Ω)`
+divisor of the classic `sin((1-t)Ω)/sin(Ω) · a + sin(tΩ)/sin(Ω) · b`
+formulation, which collapses for near- and exactly-antipodal inputs.
+Adapted from Google's S2 geometry library (see [`S2::Interpolate`](https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/s2edge_distances.cc#L77)
+and [`S2::GetPointOnLine`](https://github.com/google/s2geometry/blob/a4f0cf58a9cfc214585c39de6e3682384fac0917/src/s2/s2edge_distances.cc#L47)).
+
+For exactly antipodal `a` and `b` the great circle is mathematically
+ambiguous; `robust_cross_product` returns a deterministic perpendicular via
+its symbolic-perturbation branch, so the result is still a well-defined unit
+vector on *some* great circle through both points.
+
 ## Examples
 
 ```jldoctest
 julia> using GeometryOps.UnitSpherical
 
 julia> slerp(UnitSphericalPoint(1, 0, 0), UnitSphericalPoint(0, 1, 0), 0.5)
 3-element UnitSphericalPoint{Float64} with indices SOneTo(3):
- 0.7071067811865475
+ 0.7071067811865476
  0.7071067811865475
  0.0
 ```
 """
 function slerp(a::UnitSphericalPoint, b::UnitSphericalPoint, i01::Number)
-    if a == b
-        return a
-    end
+    i01 == 0 && return a
+    i01 == 1 && return b
+    a == b && return a
     Ω = spherical_distance(a, b)
-    sinΩ = sin(Ω)
-    return (sin((1-i01)*Ω) / sinΩ) * a + (sin(i01*Ω)/sinΩ) * b
+    dir = normalize(cross(robust_cross_product(a, b), a))
+    r = i01 * Ω
+    return normalize(cos(r) * a + sin(r) * dir)
 end
 
 function slerp(a::UnitSphericalPoint, b::UnitSphericalPoint, i01s::AbstractVector{<: Number})
-    if a == b
-        return fill(a, i01s)
-    end
+    a == b && return fill(a, size(i01s))
     Ω = spherical_distance(a, b)
-    sinΩ = sin(Ω)
-    return @. (sin((1 - i01s) * Ω) / sinΩ) * a + (sin(i01s * Ω) / sinΩ) * b
+    dir = normalize(cross(robust_cross_product(a, b), a))
+    return [begin
+        t == 0 ? a :
+        t == 1 ? b :
+        normalize(cos(t * Ω) * a + sin(t * Ω) * dir)
+    end for t in i01s]
 end
 
 #=