Spherical point-in-face and lat-lon face bounds (AccuSphGeom split from #1513)

benchmarks/geometry_kernels.py

@@ -35,6 +35,18 @@ def _unit(v):
 _X2 = _unit(np.array([0.0, 1.0, 0.3]))
 _CONST_Z = 0.3
 
+# Polygon for point-in-polygon (spherical triangle)
+_POLY = np.array(
+    [
+        _unit(np.array([1.0, 0.0, 0.1])),
+        _unit(np.array([0.0, 1.0, 0.1])),
+        _unit(np.array([-1.0, 0.0, 0.5])),
+    ],
+    dtype=np.float64,
+)
+_Q_INSIDE = _unit(np.array([0.1, 0.3, 0.9]))
+_Q_OUTSIDE = _unit(np.array([-0.5, -0.5, -0.7]))
+
 
 class EFTPrimitives:
     """Benchmark the low-level EFT building blocks: two_sum, two_prod,
@@ -205,3 +217,21 @@ def time_try_gca_const_lat_intersection(self):
     def time_gca_const_lat_intersection(self):
         """Layer 3: dispatcher (full public API)."""
         self.gca_const_lat_intersection(self.gca_cart, _CONST_Z)
+
+
+class PointInPolygonSphere:
+    """Benchmark the spherical point-in-polygon kernel."""
+
+    def setup(self):
+        from uxarray.grid.point_in_face import _point_in_polygon_sphere
+
+        self._point_in_polygon_sphere = _point_in_polygon_sphere
+
+        _point_in_polygon_sphere(_Q_INSIDE, _POLY)
+        _point_in_polygon_sphere(_Q_OUTSIDE, _POLY)
+
+    def time_point_inside(self):
+        self._point_in_polygon_sphere(_Q_INSIDE, _POLY)
+
+    def time_point_outside(self):
+        self._point_in_polygon_sphere(_Q_OUTSIDE, _POLY)

docs/userguide.rst

@@ -97,6 +97,9 @@ Supplementary Guides
 
 These user guides provide additional details about specific features in UXarray.
 
+`Accurate Spherical Geometry <user-guide/spherical-geometry-accuracy.ipynb>`_
+ How UXarray uses compensated arithmetic to avoid catastrophic cancellation in cross-product and point-in-polygon operations
+
 `Working with HEALPix Grids <user-guide/healpix.ipynb>`_
  Use UXarray with HEALPix
 
@@ -131,6 +134,7 @@ These user guides provide additional details about specific features in UXarray.
    user-guide/dual-mesh.ipynb
    user-guide/structured.ipynb
    user-guide/from-points.ipynb
+   user-guide/spherical-geometry-accuracy.ipynb
    user-guide/healpix.ipynb
    user-guide/holoviz.ipynb
    user-guide/from_file.ipynb

test/grid/geometry/test_accusphgeom_baseline.py

@@ -6,6 +6,8 @@
 Specific C++ tests mirrored here:
     tests/test_gca_gca_intersection_baseline.cpp     — 31 near-tangent GCA pairs
     tests/test_gca_constlat_intersection_baseline.cpp — 200 arc/latitude cases
+    tests/test_pip_robust.cpp                         — simple spherical triangle
+    tests/test_pip_complicated.cpp                    — 12-vertex concave polygon
 
 The C++ library uses ultra-tight tolerances (3–100 ULP) backed by Shewchuk
 adaptive precision and a geogram fallback.  This Python port implements only
@@ -14,6 +16,7 @@
 
     GCA-GCA intersection:     3e-8  (C++ reference: 1e-8)
     GCA-const-lat intersection: 1e-13 (C++ reference: 3–100 ULP ≈ 7e-16–2e-14)
+    Point-in-polygon:         exact location codes (same as C++)
 """
 
 import math
@@ -23,6 +26,13 @@
 import pytest
 
 from uxarray.grid.intersections import gca_const_lat_intersection, gca_gca_intersection
+from uxarray.grid.point_in_face import (
+    _LOC_INSIDE,
+    _LOC_ON_EDGE,
+    _LOC_ON_VERTEX,
+    _LOC_OUTSIDE,
+    _point_in_polygon_sphere,
+)
 
 _DATA_DIR = os.path.join(os.path.dirname(__file__), "data", "accusphgeom")
 _GCA_GCA_CSV = os.path.join(
@@ -133,3 +143,66 @@ def test_gca_constlat_intersection_baseline(gca_constlat_rows, idx):
     dy = result[0, 1] - by
     err = math.sqrt(dx * dx + dy * dy)
     assert err < 5e-15, f"case_id={case_id}: err_xy={err:.3e} ≥ 5e-15"
+
+
+# ── Point-in-polygon: simple spherical triangle ───────────────────────────────
+# From test_pip_robust.cpp: triangle A=(1,0,0) B=(0,1,0) C=(0,0,1)
+
+_SIMPLE_POLY = np.array(
+    [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]], dtype=np.float64
+)
+
+
+def test_pip_simple_on_vertex():
+    q = np.array([1.0, 0.0, 0.0])
+    assert _point_in_polygon_sphere(q, _SIMPLE_POLY) == _LOC_ON_VERTEX
+
+
+def test_pip_simple_on_edge():
+    # Normalize([1,1,0]) — midpoint of edge AB
+    q = np.array([0.70710678118654752, 0.70710678118654752, 0.0])
+    assert _point_in_polygon_sphere(q, _SIMPLE_POLY) == _LOC_ON_EDGE
+
+
+def test_pip_simple_inside():
+    q = np.array([1.0, 1.0, 1.0])
+    q = q / np.linalg.norm(q)
+    assert _point_in_polygon_sphere(q, _SIMPLE_POLY) == _LOC_INSIDE
+
+
+# ── Point-in-polygon: complicated 12-vertex polygon ──────────────────────────
+# From test_pip_complicated.cpp (Tier 4 / no-global-id overload)
+
+_COMPLICATED_POLY = np.array(
+    [
+        [0.77114888623389370, -0.15726142646764130, 0.61692644537707060],
+        [0.45249789144681710, -0.75061357063415830, 0.48148200985709080],
+        [0.68946150885186746, -0.59933974587969335, 0.40673664307580021],
+        [0.53398361424012150, -0.82144802877974800, 0.20021147753544170],
+        [0.72547341102583852, -0.63064441484306173, 0.27563735581699919],
+        [0.90662646752004000, -0.37288916572560260, 0.19743889808393390],
+        [0.74736479846796566, -0.64967430761889954, 0.13917310096006544],
+        [0.75468084319451650, -0.65603404827296060, -0.00872653549837396],
+        [0.49138625363591330, -0.85368085756667700, -0.17253562867386300],
+        [0.86555356123625300, -0.23932615843504300, -0.43993183849315200],
+        [0.73819995144420940, -0.26096774566031860, -0.62205841157622660],
+        [0.60166139617200880, -0.05234812405382043, -0.79703402578835670],
+    ],
+    dtype=np.float64,
+)
+
+_PIP_CASES = [
+    ([0.75367527697268680, -0.65515992289232780, -0.05233595624294383], _LOC_INSIDE, "Q1 inside"),
+    ([0.92054211727315200, -0.38498585550407840, 0.06624274592780397], _LOC_INSIDE, "Q2 inside"),
+    ([0.53882393432914170, -0.82565565483991800, 0.16721694718218960], _LOC_OUTSIDE, "Q3 outside"),
+    ([0.63494819288856630, -0.65761549896072850, 0.40544130015845230], _LOC_OUTSIDE, "Q4 outside"),
+    # Q5 is exactly vertex P8 (0-indexed)
+    ([0.49138625363591330, -0.85368085756667700, -0.17253562867386300], _LOC_ON_VERTEX, "Q5 on vertex"),
+]
+
+
+@pytest.mark.parametrize("q_xyz,expected,name", _PIP_CASES)
+def test_pip_complicated(q_xyz, expected, name):
+    q = np.array(q_xyz, dtype=np.float64)
+    result = _point_in_polygon_sphere(q, _COMPLICATED_POLY)
+    assert result == expected, f"{name}: expected {expected}, got {result}"