|
| 1 | +# # Spherical Caps |
| 2 | + |
| 3 | +#= |
| 4 | +```@meta |
| 5 | +CollapsedDocStrings = true |
| 6 | +``` |
| 7 | +
|
| 8 | +```@docs; canonical=false |
| 9 | +SphericalCap |
| 10 | +circumcenter_on_unit_sphere |
| 11 | +``` |
| 12 | +
|
| 13 | +## What is SphericalCap? |
| 14 | +
|
| 15 | +A spherical cap represents a section of a unit sphere about some point, bounded by a radius. |
| 16 | +It is defined by a center point on the unit sphere and a radius (in radians). |
| 17 | +
|
| 18 | +Spherical caps are used in: |
| 19 | +- Representing circular regions on a spherical surface |
| 20 | +- Approximating and bounding spherical geometries |
| 21 | +- Spatial indexing and filtering on the unit sphere |
| 22 | +- Implementing containment, intersection, and disjoint predicates |
| 23 | +
|
| 24 | +The `SphericalCap` type offers multiple constructors to create caps from: |
| 25 | +- UnitSphericalPoint and radius |
| 26 | +- Geographic coordinates and radius |
| 27 | +- Three points on the unit sphere (circumcircle) |
| 28 | +
|
| 29 | +## Examples |
| 30 | +
|
| 31 | +```@example sphericalcap |
| 32 | +using GeometryOps |
| 33 | +using GeoInterface |
| 34 | +
|
| 35 | +# Create a spherical cap from a point and radius |
| 36 | +point = UnitSphericalPoint(1.0, 0.0, 0.0) # Point on the unit sphere |
| 37 | +cap = SphericalCap(point, 0.5) # Cap with radius 0.5 radians |
| 38 | +``` |
| 39 | +
|
| 40 | +```@example sphericalcap |
| 41 | +# Create a spherical cap from geographic coordinates |
| 42 | +lat, lon = 40.0, -74.0 # New York City (approximate) |
| 43 | +point = GeoInterface.Point(lon, lat) |
| 44 | +cap = SphericalCap(point, 0.1) # Cap with radius ~0.1 radians |
| 45 | +``` |
| 46 | +
|
| 47 | +```@example sphericalcap |
| 48 | +# Create a spherical cap from three points (circumcircle) |
| 49 | +p1 = UnitSphericalPoint(1.0, 0.0, 0.0) |
| 50 | +p2 = UnitSphericalPoint(0.0, 1.0, 0.0) |
| 51 | +p3 = UnitSphericalPoint(0.0, 0.0, 1.0) |
| 52 | +cap = SphericalCap(p1, p2, p3) |
| 53 | +``` |
| 54 | +
|
| 55 | +=# |
| 56 | + |
| 57 | +# Spherical cap implementation |
| 58 | +struct SphericalCap{T} |
| 59 | + point::UnitSphericalPoint{T} |
| 60 | + radius::T |
| 61 | + # cartesian_radius::T # TODO: compute using cosine(radius) |
| 62 | +end |
| 63 | + |
| 64 | +SphericalCap(point::UnitSphericalPoint{T}, radius::Number) where T = SphericalCap{T}(point, convert(T, radius)) |
| 65 | +SphericalCap(point, radius::Number) = SphericalCap(GI.trait(point), point, radius) |
| 66 | + |
| 67 | +SphericalCap(geom) = SphericalCap(GI.trait(geom), geom) |
| 68 | +SphericalCap(t::GI.AbstractGeometryTrait, geom) = SphericalCap(t, geom, 0) |
| 69 | + |
| 70 | +function SphericalCap(::GI.PointTrait, point, radius::Number) |
| 71 | + return SphericalCap(UnitSphereFromGeographic()(point), radius) |
| 72 | +end |
| 73 | +# TODO: add implementations for line string and polygon traits |
| 74 | +# That will require a minimum bounding circle implementation. |
| 75 | +# TODO: add implementations for multitraits based on this |
| 76 | + |
| 77 | +# TODO: this returns an approximately antipodal point... |
| 78 | + |
| 79 | +# TODO: exact-predicate intersection |
| 80 | +# This is all inexact and thus subject to floating point error |
| 81 | +function _intersects(x::SphericalCap, y::SphericalCap) |
| 82 | + spherical_distance(x.point, y.point) <= x.radius + y.radius |
| 83 | +end |
| 84 | + |
| 85 | +_disjoint(x::SphericalCap, y::SphericalCap) = !_intersects(x, y) |
| 86 | + |
| 87 | +function _contains(big::SphericalCap, small::SphericalCap) |
| 88 | + dist = spherical_distance(big.point, small.point) |
| 89 | + # small circle fits in big circle |
| 90 | + return dist + small.radius < big.radius |
| 91 | +end |
| 92 | +function _contains(cap::SphericalCap, point::UnitSphericalPoint) |
| 93 | + spherical_distance(cap.point, point) <= cap.radius |
| 94 | +end |
| 95 | + |
| 96 | +#Comment by asinghvi: this could be transformed to GO.union |
| 97 | +function _merge(x::SphericalCap, y::SphericalCap) |
| 98 | + |
| 99 | + d = spherical_distance(x.point, y.point) |
| 100 | + newradius = (x.radius + y.radius + d) / 2 |
| 101 | + if newradius < x.radius |
| 102 | + #x contains y |
| 103 | + x |
| 104 | + elseif newradius < y.radius |
| 105 | + #y contains x |
| 106 | + y |
| 107 | + else |
| 108 | + excenter = 0.5 * (1 - (x.radius - y.radius) / d) |
| 109 | + newcenter = slerp(x.point, y.point, excenter) |
| 110 | + SphericalCap(newcenter, newradius) |
| 111 | + end |
| 112 | +end |
| 113 | + |
| 114 | +function circumcenter_on_unit_sphere(a::UnitSphericalPoint, b::UnitSphericalPoint, c::UnitSphericalPoint) |
| 115 | + LinearAlgebra.normalize( |
| 116 | + LinearAlgebra.cross(a, b) + |
| 117 | + LinearAlgebra.cross(b, c) + |
| 118 | + LinearAlgebra.cross(c, a) |
| 119 | + ) |
| 120 | +end |
| 121 | + |
| 122 | +"Get the circumcenter of the triangle (a, b, c) on the unit sphere. Returns a normalized 3-vector." |
| 123 | +function SphericalCap(a::UnitSphericalPoint, b::UnitSphericalPoint, c::UnitSphericalPoint) |
| 124 | + circumcenter = circumcenter_on_unit_sphere(a, b, c) |
| 125 | + circumradius = spherical_distance(a, circumcenter) |
| 126 | + return SphericalCap(circumcenter, circumradius) |
| 127 | +end |
| 128 | + |
| 129 | +function _is_ccw_unit_sphere(v_0::S, v_c::S, v_i::S) where S <: UnitSphericalPoint |
| 130 | + # checks if the smaller interior angle for the great circles connecting u-v and v-w is CCW |
| 131 | + return(LinearAlgebra.dot(LinearAlgebra.cross(v_c - v_0, v_i - v_c), v_i) < 0) |
| 132 | +end |
| 133 | + |
| 134 | +function angle_between(a::S, b::S, c::S) where S <: UnitSphericalPoint |
| 135 | + ab = b - a |
| 136 | + bc = c - b |
| 137 | + norm_dot = (ab ⋅ bc) / (LinearAlgebra.norm(ab) * LinearAlgebra.norm(bc)) |
| 138 | + angle = acos(clamp(norm_dot, -1.0, 1.0)) |
| 139 | + if _is_ccw_unit_sphere(a, b, c) |
| 140 | + return angle |
| 141 | + else |
| 142 | + return 2π - angle |
| 143 | + end |
| 144 | +end |
0 commit comments