Cast lasso rays in the core instead of once per polygon edge

python/xy/_native.py

@@ -24,7 +24,7 @@
 
 from .config import MAX_CONTOUR_WORK, MAX_SCREEN_DIM
 
-ABI_VERSION = 40
+ABI_VERSION = 41
 
 # Rust reports invalid arguments (and, via the ffi_guard panic shield, any
 # internal panic) by returning `usize::MAX` from size-returning entry points.
@@ -546,6 +546,18 @@ def _load() -> ctypes.CDLL:
         ctypes.c_double,
         ctypes.c_void_p,
     ]
+    lib.xy_polygon_select.restype = ctypes.c_size_t
+    lib.xy_polygon_select.argtypes = [
+        ctypes.c_void_p,  # x
+        ctypes.c_void_p,  # y
+        ctypes.c_size_t,  # len
+        ctypes.c_void_p,  # candidate row ids
+        ctypes.c_size_t,  # candidate count
+        ctypes.c_void_p,  # polygon x
+        ctypes.c_void_p,  # polygon y
+        ctypes.c_size_t,  # polygon vertices
+        ctypes.c_void_p,  # output row IDs
+    ]
     lib.xy_sample_mask.restype = ctypes.c_int32
     lib.xy_sample_mask.argtypes = [
         ctypes.c_void_p,
@@ -2510,6 +2522,47 @@ def range_indices(
     return out if written == len(out) else out[:written].copy()
 
 
+def polygon_select(
+    x: npt.NDArray[np.float64],
+    y: npt.NDArray[np.float64],
+    rows: npt.NDArray[np.uint32],
+    poly_x: npt.NDArray[np.float64],
+    poly_y: npt.NDArray[np.float64],
+) -> npt.NDArray[np.uint32]:
+    """Those of `rows` whose canonical point is inside the lasso polygon.
+
+    Even-odd ray casting, order preserved. Non-finite coordinates are never
+    inside, matching `range_indices`."""
+    x = _as_f64(x, "x")
+    y = _as_f64(y, "y")
+    if len(x) != len(y):
+        raise ValueError("x and y must have equal length")
+    poly_x = _as_f64(poly_x, "polygon x")
+    poly_y = _as_f64(poly_y, "polygon y")
+    if len(poly_x) != len(poly_y):
+        raise ValueError("polygon x and y must have equal length")
+    rows = np.ascontiguousarray(rows, dtype=np.uint32)
+    if rows.ndim != 1:
+        raise ValueError("rows must be 1-D")
+    out = np.empty(len(rows), dtype=np.uint32)
+    if len(rows) == 0:
+        return out
+    written = _lib.xy_polygon_select(
+        _ptr_f64(x),
+        _ptr_f64(y),
+        len(x),
+        rows.ctypes.data,
+        len(rows),
+        _ptr_f64(poly_x),
+        _ptr_f64(poly_y),
+        len(poly_x),
+        out.ctypes.data,
+    )
+    if written == _USIZE_MAX:
+        raise ValueError("invalid polygon_select arguments")
+    return out if written == len(out) else out[:written].copy()
+
+
 def sample_mask(
     ids: npt.NDArray[np.uint64],
     seed: int,

python/xy/interaction.py

@@ -225,6 +225,14 @@ def select_polygon(
 
     The polygon's bounding box first reuses the zone-pruned range predicate;
     ray casting then runs only on those candidates rather than every row.
+
+    The cast itself is native (`kernels.polygon_select`). Vectorizing it over
+    NumPy meant one pass per polygon EDGE, each materializing full-length
+    crossing and edge-intersection temporaries: a 64-vertex lasso over 160k
+    candidates moved roughly half a gigabyte through memory to answer a
+    question that fits in registers. Crossing parity is order-independent, so
+    the native kernel walks edges inside the point loop and returns the same
+    rows.
     """
     polygon = np.asarray(points, dtype=np.float64)
     if polygon.ndim != 2 or polygon.shape[1:] != (2,) or not 3 <= len(polygon) <= 2048:
@@ -239,25 +247,15 @@ def select_polygon(
         float(polygon[:, 1].max()),
         trace_id,
     )
+    poly_x = np.ascontiguousarray(polygon[:, 0])
+    poly_y = np.ascontiguousarray(polygon[:, 1])
     out: dict[int, np.ndarray] = {}
     for tid, rows in candidates.items():
         if len(rows) == 0:
             out[tid] = rows
             continue
         trace = fig.traces[tid]
-        x = trace.x.values[rows]
-        y = trace.y.values[rows]
-        inside = np.zeros(len(rows), dtype=np.bool_)
-        j = len(polygon) - 1
-        for i in range(len(polygon)):
-            xi, yi = polygon[i]
-            xj, yj = polygon[j]
-            crosses = (yi > y) != (yj > y)
-            with np.errstate(divide="ignore", invalid="ignore"):
-                edge_x = (xj - xi) * (y - yi) / (yj - yi) + xi
-            inside ^= crosses & (x < edge_x)
-            j = i
-        out[tid] = rows[inside]
+        out[tid] = kernels.polygon_select(trace.x.values, trace.y.values, rows, poly_x, poly_y)
     return out
 
 

python/xy/kernels.py

@@ -62,6 +62,7 @@
 valid_indices_f64 = _impl.valid_indices_f64
 remap_u8 = _impl.remap_u8
 range_indices = _impl.range_indices
+polygon_select = _impl.polygon_select
 sample_mask = _impl.sample_mask
 sample_range_indices = _impl.sample_range_indices
 stratified_sample_range_u8 = _impl.stratified_sample_range_u8
@@ -121,6 +122,7 @@
     "marching_triangles",
     "min_max",
     "normalize_f32",
+    "polygon_select",
     "polygon_triangles",
     "pyramid_append",
     "pyramid_build",

spec/design/rust-engine.md

@@ -25,7 +25,7 @@ clear ImportError when it can't load, with no pure-Python fallback.
 
 | Concern | Today | Verdict |
 |---|---|---|
-| zone maps, encode_f32, m4, bin_2d, bin_2d_mean_color, min_max, histogram_uniform, normalize_f32, range/validity indices, local_log_density | Rust (ABI v40) | correct — new equal-length x/y columns use a paired zone-map call with bit-identical per-column reductions; full-domain density first paint fuses binning with uniform or counted-u8 overlay sampling while retaining exact standalone outputs; mean-color binning (LOD doc §2) is an integer-only pipeline (checked-in sRGB⇄linear-u16 tables, alpha-weighted u64 sums) so grids are bitwise deterministic across thread counts and platforms; mesh/rectangle validity scans consume only columns not already proven finite by zone metadata |
+| zone maps, encode_f32, m4, bin_2d, bin_2d_mean_color, min_max, histogram_uniform, normalize_f32, range/validity indices, polygon (lasso) selection, local_log_density | Rust (ABI v41) | correct — new equal-length x/y columns use a paired zone-map call with bit-identical per-column reductions; lasso ray casting (`xy_polygon_select`, §34) walks edges inside the point loop instead of one NumPy pass per edge, and buckets edges by y so a point tests only those spanning its row — the answer is the same crossing parity, without the per-edge full-length temporaries that made a 2048-vertex lasso cost ~370 ms over 160k candidates. Bucketing is declined below 16 vertices (the build outweighs the scans it saves) and for a non-finite polygon (slab bounds stop meaning anything); above 16 the slab count tracks the vertex count, capped at 256 and shed further so the CSR never exceeds 2²⁰ entries, which bounds both the allocation and its u32 cursors independently of what a caller passes; full-domain density first paint fuses binning with uniform or counted-u8 overlay sampling while retaining exact standalone outputs; mean-color binning (LOD doc §2) is an integer-only pipeline (checked-in sRGB⇄linear-u16 tables, alpha-weighted u64 sums) so grids are bitwise deterministic across thread counts and platforms; mesh/rectangle validity scans consume only columns not already proven finite by zone metadata |
 | fixed-width string/bytes/bool factorization | Rust (ABI v36) | correct — compact palettes use a bounded L1-resident codebook with full-record collision checks and emit exact counts; U1 uses a direct Unicode-scalar table with endian support; ≥512k rows probe a prefix then encode disjoint chunks in parallel, merging late labels by canonical first-row order before any retry; Python sees only unique labels and retains display-label ordering policy |
 | static display-list raster, row-banded polyline/point/segment paint, batched fill+stroke triangle meshes, affine scatter projection plus typed color/size resolution, density/heatmap colormap and sampling | Rust (ABI v36) | correct — commands borrow f32/u8 payload or canonical spans synchronously; compact stratified sampling reuses factorization counts; batched/banded output is byte-identical |
 | signal processing: `xy_rfft`, `xy_welch_spectra`, `xy_spectrogram` | Rust (ABI v36) | correct — O(N) transforms over sample columns; window/segment policy stays in Python |