"""Lines spatial ordering (dual: vertices in D-space + segments in 2D-space).
Line vertices are ordered like Points (in D-space); segments are ordered in
(2xD)-space by concatenating both endpoints, so spatially-coherent segments land
in the same chunk. Chunk bounds include the line-width extent on spatial axes.
"""
from __future__ import annotations
from typing import Literal, Optional
import numpy as np
from numpy.typing import NDArray
from luxar.core import Dimension
from .bounds import (
_BARRIER_BOUND_EPS,
_normalise_coord_slack,
_normalise_scalar_slack,
_normalise_slice_dims,
_store_outward_f32_array,
)
from .curves.hilbert import hilbert_encode_nd
from .curves.morton import morton_encode_128bit, morton_encode_nd
from .grid import normalize_coords_to_grid
from .points import sort_points_compound
[docs]
def convert_to_indexed(
n_vertices: int, line_type: str, indices: Optional[np.ndarray]
) -> np.ndarray:
"""Convert any line type to indexed segment pairs.
All line types are internally converted to the unified indexed representation
for efficient spatial ordering and storage.
Args:
n_vertices: Number of vertices in the lines
line_type: One of "segments", "polyline", "loop", "indexed"
indices: For "indexed" type, the user-provided index array
Returns:
segments: (S, 2) uint32 array of vertex index pairs
Raises:
ValueError: If line_type is invalid or indices missing for "indexed" type
"""
if line_type == "indexed":
if indices is None:
raise ValueError("Indexed line type requires indices array")
return indices.reshape(-1, 2).astype(np.uint32)
elif line_type == "segments":
# Each pair of consecutive vertices forms a segment
return np.arange(n_vertices, dtype=np.uint32).reshape(-1, 2)
elif line_type == "polyline":
# Connect consecutive vertices: (0,1), (1,2), (2,3), ...
return np.column_stack(
[
np.arange(n_vertices - 1, dtype=np.uint32),
np.arange(1, n_vertices, dtype=np.uint32),
]
)
elif line_type == "loop":
# Like polyline but connect last to first
return np.column_stack(
[
np.arange(n_vertices, dtype=np.uint32),
np.roll(np.arange(n_vertices, dtype=np.uint32), -1),
]
)
else:
raise ValueError(
f"Invalid line_type '{line_type}'. Must be one of: segments, polyline, loop, indexed"
)
[docs]
def sort_segments_compound(
segment_coords_2d: np.ndarray,
dimensions: list[Dimension],
method: Literal["morton", "hilbert"] = "hilbert",
) -> tuple[np.ndarray, dict]:
"""Sort segments using compound ordering in (2×D)-space.
Segments are represented as (2×D)-dimensional points by concatenating both
endpoint coordinates. This captures the full geometric nature of segments
(position, orientation, length) for spatial coherence.
The compound ordering strategy is:
- Primary sort: Discrete dimensions from both endpoints (lexicographic)
- Secondary sort: Morton/Hilbert code of spatial dimensions from both endpoints
Args:
segment_coords_2d: Segment coordinates in (2×D)-space, shape (S, 2*d)
First d columns are start point, second d columns are end point.
dimensions: Original Dimension objects (may have more dims than data)
method: Spatial curve method ("morton" or "hilbert"), default "hilbert"
Returns:
sort_indices: Indices to reorder segments
metadata: Dict with ordering metadata in (2×D)-space
"""
n_segments, n_dims_2d = segment_coords_2d.shape
# Use actual data dimensionality (segment_coords_2d is 2×D, so D = n_dims_2d // 2)
n_dims_original = n_dims_2d // 2
# Use only the first n_dims_original dimensions from the dimensions list
dimensions = dimensions[:n_dims_original]
# Build (2×D) dimension classification
# Discrete dims from both endpoints: if dim 3,4 are discrete in D-space,
# then dims 3,4,d+3,d+4 are discrete in (2×D)-space
slice_dims_2d = []
ordering_dims_2d = []
for i, d in enumerate(dimensions):
if d.discrete and not d.display:
slice_dims_2d.append(i) # Start point discrete dim
slice_dims_2d.append(n_dims_original + i) # End point discrete dim
else:
ordering_dims_2d.append(i) # Start point spatial dim
ordering_dims_2d.append(n_dims_original + i) # End point spatial dim
# Compute bits per ordering dimension (auto-select 128-bit if needed)
if ordering_dims_2d:
bits_per_dim = min(21, 64 // len(ordering_dims_2d))
use_128bit = bits_per_dim < 10
if use_128bit:
bits_per_dim = min(21, 128 // len(ordering_dims_2d))
else:
bits_per_dim = 21
use_128bit = False
# Extract ordering dimension coordinates
spatial_codes: tuple[np.ndarray, np.ndarray] | np.ndarray
if ordering_dims_2d:
ordering_coords = segment_coords_2d[:, ordering_dims_2d]
ordering_min = ordering_coords.min(axis=0)
ordering_max = ordering_coords.max(axis=0)
# Normalize to grid
grid_coords = normalize_coords_to_grid(
ordering_coords, ordering_min, ordering_max, 2**bits_per_dim
)
# Compute spatial curve codes
if method == "morton":
if use_128bit:
high, low = morton_encode_128bit(grid_coords, bits_per_dim)
spatial_codes = (high, low) # Tuple for lexsort
else:
spatial_codes = morton_encode_nd(grid_coords, bits_per_dim)
elif method == "hilbert":
# Hilbert doesn't have 128-bit implementation - fall back to morton for high dims
if use_128bit:
high, low = morton_encode_128bit(grid_coords, bits_per_dim)
spatial_codes = (high, low)
else:
spatial_codes = hilbert_encode_nd(grid_coords, bits_per_dim)
else:
raise ValueError(f"Unknown method: {method}")
else:
# No ordering dimensions - all discrete
spatial_codes = np.zeros(n_segments, dtype=np.uint64)
ordering_min = np.array([])
ordering_max = np.array([])
# Create compound sort key
if slice_dims_2d:
# Extract discrete dimension values
slice_values = segment_coords_2d[:, slice_dims_2d]
# Create sort keys: (discrete_tuple, spatial_code)
if isinstance(spatial_codes, tuple):
# 128-bit: lexsort by (low, high, discrete_dims...)
high, low = spatial_codes
sort_indices = np.lexsort(
[low, high]
+ [slice_values[:, i] for i in range(len(slice_dims_2d) - 1, -1, -1)]
)
else:
sort_indices = np.lexsort(
[spatial_codes]
+ [slice_values[:, i] for i in range(len(slice_dims_2d) - 1, -1, -1)]
)
else:
# Pure spatial ordering (no discrete dimensions)
if isinstance(spatial_codes, tuple):
high, low = spatial_codes
sort_indices = np.lexsort([low, high])
else:
sort_indices = np.argsort(spatial_codes)
# Build metadata
metadata = {
"ordering": method,
"slice_dims": slice_dims_2d,
"ordering_dims": ordering_dims_2d,
"ordering_min": ordering_min.tolist() if len(ordering_min) > 0 else [],
"ordering_max": ordering_max.tolist() if len(ordering_max) > 0 else [],
"ordering_bits_per_dim": bits_per_dim,
}
return sort_indices, metadata
[docs]
def order_lines_spatial(
vertices: np.ndarray,
segments: np.ndarray,
dimensions: list[Dimension],
method: Literal["morton", "hilbert"] = "hilbert",
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, dict]:
"""Apply dual spatial ordering to lines using Morton or Hilbert curves.
Lines have dual ordering:
1. Vertices are ordered in D-space (like Points)
2. Segments are ordered in (2×D)-space (concatenating both endpoints)
This is the main entry point for Lines spatial ordering.
Args:
vertices: (V, D) float32 vertex positions
segments: (S, 2) uint32 index pairs (from convert_to_indexed)
dimensions: List of Dimension objects (may have more dims than data)
method: Ordering method ("morton" or "hilbert")
Returns:
sorted_vertices: (V, D) float32 reordered vertex positions
sorted_segments: (S, 2) uint32 reordered and remapped segment indices
vertex_sort_indices: Indices to recover original vertex order
segment_sort_indices: Indices to recover original segment order
metadata: Dict with "vertex_ordering" and "segment_ordering" sub-dicts
"""
_V, D = vertices.shape # V unused but kept for reference
# Use only the first D dimensions from the dimensions list
# This handles cases where scene has more dimensions than data
dimensions = dimensions[:D]
# 1. Order vertices in D-space (with compound ordering for discrete dims)
vertex_sort_indices, vertex_metadata = sort_points_compound(
vertices, dimensions, method=method
)
sorted_vertices = vertices[vertex_sort_indices]
# 2. Create inverse mapping for index remapping
# After sorting, vertex at original index i is now at inverse_map[i]
inverse_map = np.argsort(vertex_sort_indices).astype(np.uint32)
remapped_segments = inverse_map[segments]
# 3. Build (2×D) segment coordinates for ordering
# Each segment becomes a 2D point: [x1, y1, z1, x2, y2, z2] for 3D
segment_coords_2d = np.concatenate(
[
sorted_vertices[remapped_segments[:, 0]], # Start points
sorted_vertices[remapped_segments[:, 1]], # End points
],
axis=1,
) # Shape: (S, 2×D)
# 4. Order segments in (2×D)-space (with compound ordering)
segment_sort_indices, segment_metadata = sort_segments_compound(
segment_coords_2d, dimensions, method=method
)
sorted_segments = remapped_segments[segment_sort_indices]
return (
sorted_vertices,
sorted_segments,
vertex_sort_indices,
segment_sort_indices,
{"vertex_ordering": vertex_metadata, "segment_ordering": segment_metadata},
)
[docs]
def compute_vertex_chunk_bounds(
vertices: np.ndarray,
chunk_size: int,
slice_dims: Optional[list[int]] = None,
*,
coord_slack: Optional[NDArray[np.float64]] = None,
) -> np.ndarray:
"""Compute chunk bounding boxes for vertices (no radius/width expansion).
Spatial axes carry no pad (a vertex is a point), so the interval is exact;
the barrier epsilon, however, is a small ABSOLUTE quantity. Both are computed
in float64 and narrowed to the float32 store with OUTWARD rounding (see
:func:`_store_outward_f32_array`), so a stored bound is never tighter than
the footprint at any coordinate magnitude.
QUANTISATION SLACK (``coord_slack``): that holds for the vertices AS HANDED
IN. Under the default AUTO encoding the vertices are stored as per-axis
uint16 fixed point, so a DECODED vertex can sit up to half a quantum
(``extent/131070``) outside a bound derived from the authored one on a
non-gridded axis. ``coord_slack`` is that displacement, per axis, added
outward on EVERY dimension — including on top of ``_BARRIER_BOUND_EPS`` on a
barrier axis. The compiler supplies it from the encoder's own predicate; a
direct caller that omits it gets bounds for the authored vertices only. See
:func:`~luxar.io._ordering.points.compute_chunk_bounds_points`.
Args:
vertices: Vertex positions (already sorted), shape (V, D)
chunk_size: Number of vertices per chunk
slice_dims: Indices of discrete (non-spatial) dimensions (padded by a
float-boundary epsilon only — the reader's per-dimension tolerance
owns the query reach). An index outside ``[0, D)`` raises
``ValueError`` (see :func:`_normalise_slice_dims`).
coord_slack: Per-axis outward pad, shape ``(D,)``, covering how far the
store can move a vertex from the value passed here. Default None ⇒
zero on every axis. Validated by :func:`_normalise_coord_slack`.
Returns:
chunk_bounds: Bounding boxes, shape (num_chunks, D, 2)
"""
n_vertices, n_dims = vertices.shape
num_chunks = (n_vertices + chunk_size - 1) // chunk_size
discrete_dims = _normalise_slice_dims(slice_dims, n_dims)
slack = _normalise_coord_slack(coord_slack, n_dims)
chunk_bounds = np.zeros((num_chunks, n_dims, 2), dtype=np.float32)
for chunk_idx in range(num_chunks):
start_idx = chunk_idx * chunk_size
end_idx = min(start_idx + chunk_size, n_vertices)
chunk_verts = vertices[start_idx:end_idx]
# Reduce per DIMENSION, not with a single ``min(axis=0)``: numpy's
# outer-axis reduce over a 3-or-4-element inner row is several times
# slower than one reduce per column, and these run over every chunk of
# every dataset. Only the outward STORE is vectorised.
mins = np.empty(n_dims, dtype=np.float64)
maxs = np.empty(n_dims, dtype=np.float64)
for d in range(n_dims):
col = chunk_verts[:, d]
mins[d] = col.min()
maxs[d] = col.max()
for d in discrete_dims:
# Discrete: tight bounds padded by a float-boundary epsilon only
# (the reader's per-dimension tolerance owns the query reach).
mins[d] -= _BARRIER_BOUND_EPS
maxs[d] += _BARRIER_BOUND_EPS
# Spatial dims keep their exact min/max; the outward store only ever
# moves a bound by the one ULP the cast itself would have swallowed.
# The quantisation displacement is added first, in float64: the outward
# store cannot recover a pad float32 arithmetic already rounded away.
lo32, hi32 = _store_outward_f32_array(mins - slack, maxs + slack)
chunk_bounds[chunk_idx, :, 0] = lo32
chunk_bounds[chunk_idx, :, 1] = hi32
return chunk_bounds
[docs]
def compute_segment_chunk_bounds(
vertices: np.ndarray,
segments: np.ndarray,
widths: np.ndarray,
chunk_size: int,
slice_dims: Optional[list[int]] = None,
*,
coord_slack: Optional[NDArray[np.float64]] = None,
scalar_slack: Optional[float] = None,
) -> np.ndarray:
"""Compute chunk bounding boxes for segments (includes line width).
Segment bounds are in D-dimensional space (not 2×D) for view frustum
intersection tests. Each segment's bounds include the line width extent.
The width interval is accumulated in float64 and narrowed to the float32
store with OUTWARD rounding (see :func:`_store_outward_f32_array`), so a
stored bound is never tighter than the footprint of the AUTHORED widths at
any coordinate magnitude — not only where a small width happens to survive
float32 arithmetic and a round-to-nearest store.
SCALAR QUANTISATION SLACK (``scalar_slack``): ``widths`` is itself a
POSITIVE_SCALAR and may decode larger than the authored value. This single
per-array pad is added to the width on SPATIAL dimensions only; barrier
dimensions still receive no footprint expansion. The compiler supplies it
from
:meth:`~luxar.encoding.encoder.ArrayEncoder.positive_scalar_round_trip_slack`;
a direct caller that omits it gets authored-width bounds.
QUANTISATION SLACK (``coord_slack``): that holds for the vertices AS HANDED
IN. Under the default AUTO encoding the vertices are stored as per-axis
uint16 fixed point, so a DECODED endpoint can sit up to half a quantum
(``extent/131070``) outside a bound derived from the authored one on a
non-gridded axis. ``coord_slack`` is that displacement, per axis, added
outward on EVERY dimension — on top of the width on a spatial axis and on
top of ``_BARRIER_BOUND_EPS`` on a barrier axis. The compiler supplies it
from the encoder's own predicate (the SAME vector it gives
:func:`compute_vertex_chunk_bounds`: both builders are in D-space over the
same vertices array); a direct caller that omits it gets bounds for the
authored vertices only. See
:func:`~luxar.io._ordering.points.compute_chunk_bounds_points`.
IMPORTANT: widths must be a full (V,) array. Broadcast widths should be
expanded with np.full(V, width_value) before calling this function.
Args:
vertices: Vertex positions (already sorted), shape (V, D)
segments: Segment index pairs (already sorted), shape (S, 2)
widths: Vertex widths (already sorted), shape (V,)
chunk_size: Number of segments per chunk
slice_dims: Indices of discrete (non-spatial) dimensions (padded by a
float-boundary epsilon only — the reader's per-dimension tolerance
owns the query reach). An index outside ``[0, D)`` raises
``ValueError`` (see :func:`_normalise_slice_dims`).
coord_slack: Per-axis outward pad, shape ``(D,)``, covering how far the
store can move a vertex from the value passed here. Default None ⇒
zero on every axis. Validated by :func:`_normalise_coord_slack`.
scalar_slack: Outward pad covering how far the stored width can exceed
the authored width. Applied on spatial dimensions only; overflow
produces conservative infinite spatial bounds.
Returns:
chunk_bounds: Bounding boxes, shape (num_chunks, D, 2)
"""
_V, D = vertices.shape
S = segments.shape[0]
num_chunks = (S + chunk_size - 1) // chunk_size
discrete_dims = _normalise_slice_dims(slice_dims, D)
slack = _normalise_coord_slack(coord_slack, D)
footprint_slack = _normalise_scalar_slack(scalar_slack)
chunk_bounds = np.zeros((num_chunks, D, 2), dtype=np.float32)
for chunk_idx in range(num_chunks):
start_idx = chunk_idx * chunk_size
end_idx = min(start_idx + chunk_size, S)
chunk_segs = segments[start_idx:end_idx]
# Get vertex positions and widths for this chunk's segments. The width
# is widened to float64 so every ``coord ± width`` below is a float64
# sum: a small width added to a large coordinate is lost outright in
# float32 arithmetic, before the outward store can rescue it.
p1 = vertices[chunk_segs[:, 0]]
p2 = vertices[chunk_segs[:, 1]]
w1 = widths[chunk_segs[:, 0]].astype(np.float64, copy=False)
w2 = widths[chunk_segs[:, 1]].astype(np.float64, copy=False)
# Overflow to +inf only widens the bound, so it is conservative.
with np.errstate(over="ignore"):
max_w = np.maximum(w1, w2) + footprint_slack
# Reduce per DIMENSION, not with a single ``min(axis=0)``: numpy's
# outer-axis reduce over a 3-or-4-element inner row is several times
# slower than one reduce per column, and these run over every chunk of
# every dataset. Only the outward STORE below is vectorised.
mins = np.empty(D, dtype=np.float64)
maxs = np.empty(D, dtype=np.float64)
for d in range(D):
c1 = p1[:, d]
c2 = p2[:, d]
if d in discrete_dims:
# Discrete: no width expansion, only a float-boundary epsilon.
# The reader's per-dimension tolerance owns the query reach
# (see _BARRIER_BOUND_EPS). float() keeps the epsilon in
# float64 — a float32 column min would demote the subtraction.
mins[d] = min(float(c1.min()), float(c2.min())) - _BARRIER_BOUND_EPS
maxs[d] = max(float(c1.max()), float(c2.max())) + _BARRIER_BOUND_EPS
else:
# Spatial: include the full endpoint width extent (no /2).
mins[d] = min((c1 - max_w).min(), (c2 - max_w).min())
maxs[d] = max((c1 + max_w).max(), (c2 + max_w).max())
# Quantisation displacement on every dim, in float64, before the store.
lo32, hi32 = _store_outward_f32_array(mins - slack, maxs + slack)
chunk_bounds[chunk_idx, :, 0] = lo32
chunk_bounds[chunk_idx, :, 1] = hi32
return chunk_bounds