Skip to content

experimental UnitSpherical module that can express geometry on unit sphere #285

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 17 commits into
base: main
Choose a base branch
from

Conversation

asinghvi17
Copy link
Member

@asinghvi17 asinghvi17 commented Apr 3, 2025

TL;DR

Wrote a new UnitSpherical.jl module (experimental) for working with points on the unit sphere and spherical caps.

What changed?

Created a new file UnitSpherical.jl that implements:

  • UnitSphericalPoint struct for representing points on the unit sphere (S²) using Cartesian coordinates in ℝ³
  • SphericalCap struct for representing spherical caps
  • Coordinate transformations between geographic (lat/long) and unit spherical coordinates
  • GeoInterface compatibility for the new types
  • Utility functions for spherical geometry:
    • spherical_distance for calculating distances between points
    • slerp for spherical linear interpolation
    • circumcenter_on_unit_sphere for finding the center of a circle passing through three points: still broken
    • angle_between for calculating angles between points
    • randsphere for generating random points on the unit sphere

Quick start

  1. Import the module and create unit spherical points:

    using GeometryOps.UnitSpherical
    
    # Create a point from Cartesian coordinates
    p1 = UnitSphericalPoint(1.0, 0.0, 0.0)
    
    # Create a point from geographic coordinates (lat/long)
    p2 = UnitSphericalPoint([45.0, 90.0])
    
    # Calculate spherical distance
    dist = spherical_distance(p1, p2)
  2. Test spherical caps and operations:

    cap1 = SphericalCap(p1, 0.5)
    cap2 = SphericalCap(p2, 0.3)
    
    # Test intersection
    _intersects(cap1, cap2)
  3. Generate random points on the sphere:

    points = rand(UnitSphericalPoint, 100)

This should provide a nice foundation to start looking at R-trees that are backed by spherical caps, not extents, in preparation for @meggart's work involving this. With the new changes to as/trees SpatialTreeInterface is now representation agnostic, meaning that as long as the predicate you passed in can consume the extent like object you have chosen, any extent like thing (accessible by SpatialTreeInterface.node_extent) can be used by the STI algorithms.

Copy link
Member Author

asinghvi17 commented Apr 3, 2025

This stack of pull requests is managed by Graphite. Learn more about stacking.

@asinghvi17 asinghvi17 added enhancement New feature or request low priority Low priority issue, someone can get to it whenever they feel like it! labels Apr 3, 2025 — with Graphite App
@asinghvi17 asinghvi17 marked this pull request as ready for review April 3, 2025 20:15
@asinghvi17 asinghvi17 force-pushed the as/spherical-cap-str branch from 3d14dc5 to 6118f5e Compare April 3, 2025 20:19
@asinghvi17 asinghvi17 changed the base branch from as/trees to graphite-base/285 April 4, 2025 02:24
@asinghvi17 asinghvi17 force-pushed the as/spherical-cap-str branch from 16cf9b4 to 00ff75c Compare April 4, 2025 02:26
@asinghvi17 asinghvi17 changed the base branch from graphite-base/285 to as/trees April 4, 2025 02:26
@asinghvi17
Copy link
Member Author

For this to work well, we'll probably want to port a lot of the stuff from s2, like https://github.com/google/s2geometry/blob/58de4ea1e2f8a294e0c072c602c22232fd1433ad/src/s2/s1chord_angle.h - it's Apache 2.0 licensed so we can directly port it.

They have a good framework to handle this stuff, so if we can get the basics sorted then we're mostly there I believe. That would also be a good foundation to build an s2 indexing framework in Julia, since we would have all the basic components already.

@asinghvi17 asinghvi17 changed the base branch from as/trees to graphite-base/285 April 16, 2025 20:38
@asinghvi17 asinghvi17 force-pushed the graphite-base/285 branch from f33c8e4 to 2ca2a3c Compare May 4, 2025 16:07
@asinghvi17 asinghvi17 force-pushed the as/spherical-cap-str branch from 747b248 to 6e8baa1 Compare May 4, 2025 16:07
@asinghvi17 asinghvi17 changed the base branch from graphite-base/285 to main May 4, 2025 16:07
asinghvi17 and others added 4 commits May 4, 2025 12:35
@asinghvi17 asinghvi17 force-pushed the as/spherical-cap-str branch from 93db6d6 to 0df0f74 Compare May 5, 2025 01:48
@asinghvi17 asinghvi17 added unitspherical and removed low priority Low priority issue, someone can get to it whenever they feel like it! labels May 9, 2025
Ported the implementation from s2.  This is the core component in spherical intersection_point which we will need for polygon intersection.  Thankfully with a few tweaks, FosterHormannClipping should permit spherical points as well and I don't think we need to make any changes after that.  Predicates should also work in 3d, maybe with an implicit (0, 0, 0) point to round out the plane.

function _is_ccw_unit_sphere(v_0::S, v_c::S, v_i::S) where S <: UnitSphericalPoint
# checks if the smaller interior angle for the great circles connecting u-v and v-w is CCW
return(LinearAlgebra.dot(LinearAlgebra.cross(v_c - v_0,v_i - v_c), v_i) < 0)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This needs to use robust_cross_product, probably.



function circumcenter_on_unit_sphere(a::UnitSphericalPoint, b::UnitSphericalPoint, c::UnitSphericalPoint)
LinearAlgebra.normalize(a × b + b × c + c × a)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So does this

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request unitspherical
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants