"""Native building-shell placement and vertical-column sampling.
These operations consume source triangles, not footprint extrusions. Z remains
in its declared source datum until an explicit common building offset is applied.
"""
from __future__ import annotations
from functools import lru_cache
import numpy as np
import shapely
from pyproj import CRS, Transformer
from openworld_radio_twin.geodesy import LocalFrame
from openworld_radio_twin.models import BuildingGeometry
# Upper bound on simultaneously evaluated (triangle, grid point) pairs.
_PAIR_CHUNK = 2_000_000
_LEVEL_TOLERANCE = 1e-6
@lru_cache(maxsize=32)
def _horizontal_transform(crs: str) -> Transformer:
return Transformer.from_crs(CRS(crs).to_2d(), "EPSG:4326", always_xy=True)
[docs]
def geographic_vertices(geometry: BuildingGeometry) -> np.ndarray:
"""Return WGS84 XY with untouched source-datum Z (not ellipsoidal altitude)."""
vertices = np.asarray(geometry.vertices, dtype=np.float64)
longitude, latitude = _horizontal_transform(geometry.crs).transform(
vertices[:, 0], vertices[:, 1]
)
return np.column_stack((longitude, latitude, vertices[:, 2]))
[docs]
def horizontal_enu_vertices(geometry: BuildingGeometry, frame: LocalFrame) -> np.ndarray:
vertices = geographic_vertices(geometry)
east, north, _ = frame.transformer.transform(
vertices[:, 0], vertices[:, 1], np.zeros(len(vertices))
)
return np.column_stack((east, north, vertices[:, 2]))
[docs]
def has_ground_faces(geometry: BuildingGeometry) -> bool:
"""Whether the source shell carries GroundSurface faces to anchor on."""
return "GroundSurface" in geometry.surface_types
[docs]
def ground_projection(vertices: np.ndarray, geometry: BuildingGeometry):
"""Return the source ground contact used for anchoring.
GroundSurface faces exclude roof overhangs. A shell without them, which 3DBAG emits
for some parts, falls back to the projection of its whole shell so the building is
still placed; the caller records that rule instead of dropping the building.
"""
triangles = np.asarray(geometry.triangles, dtype=np.intp).reshape(-1, 3)
if has_ground_faces(geometry):
kinds = np.asarray(geometry.surface_types, dtype=object)
triangles = triangles[kinds == "GroundSurface"]
polygons = projected_triangles(np.asarray(vertices, dtype=np.float64), triangles)
if polygons.size == 0:
raise ValueError("Native building geometry has no faces with a ground projection")
return shapely.union_all(polygons)
[docs]
def projected_triangles(vertices: np.ndarray, triangles: np.ndarray) -> np.ndarray:
"""Vertical projections of triangles with non-negligible area, as one geometry array."""
triangles = np.asarray(triangles, dtype=np.intp).reshape(-1, 3)
if triangles.size == 0:
return np.empty(0, dtype=object)
polygons = shapely.polygons(np.asarray(vertices, dtype=np.float64)[triangles][:, :, :2])
return polygons[shapely.area(polygons) > 1e-10]
[docs]
def column_intervals(
vertices: np.ndarray,
triangles: np.ndarray,
east_axis: np.ndarray,
north_axis: np.ndarray,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""Sample solid occupancy, including compound shells and enclosed cavities.
Independent closed components can overlap or share a ground elevation. Pair
intersections within each component before combining signed interior intervals;
globally deduplicating elevations loses the multiplicity at shared surfaces.
Outward-oriented components add solid volume; inward-oriented ones bound voids.
"""
vertices = np.asarray(vertices, dtype=np.float64)
triangles = np.asarray(triangles, dtype=np.intp)
components = _shell_components(vertices, triangles)
if len(components) <= 1:
return _shell_intervals(vertices, triangles, east_axis, north_axis)
width = len(east_axis)
keys: list[np.ndarray] = []
levels: list[np.ndarray] = []
changes: list[np.ndarray] = []
for faces in components:
rows, columns, lower, upper = _shell_intervals(vertices, faces, east_axis, north_axis)
points = vertices[faces]
points = points - points.mean(axis=(0, 1)) # Stable signed volume in local coordinates.
volume = np.einsum("ij,ij->i", points[:, 0], np.cross(points[:, 1], points[:, 2])).sum()
if abs(volume) < 1e-12:
raise ValueError("Native shell component has zero signed volume")
sign = 1 if volume > 0 else -1
component_keys = rows.astype(np.int64) * width + columns
keys.extend((component_keys, component_keys))
levels.extend((lower, upper))
changes.extend(
(np.full(rows.size, sign, dtype=np.int64), np.full(rows.size, -sign, dtype=np.int64))
)
key_array = np.concatenate(keys) if keys else np.empty(0, dtype=np.int64)
if key_array.size == 0:
return (np.empty(0, dtype=np.intp), np.empty(0, dtype=np.intp), np.empty(0), np.empty(0))
level_array = np.concatenate(levels)
change_array = np.concatenate(changes)
order = np.lexsort((change_array, level_array, key_array))
result = []
for key, start, stop in _groups(key_array[order]):
crossings = zip(
level_array[order][start:stop], change_array[order][start:stop], strict=True
)
merged: list[list[float]] = []
for z, change in crossings:
if merged and z - merged[-1][0] <= _LEVEL_TOLERANCE:
merged[-1][1] += change
else:
merged.append([float(z), int(change)])
winding, interval_start = 0, None
for z, change in merged:
next_winding = winding + change
if winding == 0 and next_winding != 0:
interval_start = z
elif winding != 0 and next_winding == 0 and z - interval_start > _LEVEL_TOLERANCE:
result.append((key // width, key % width, interval_start, z))
winding = next_winding
if not result:
return (np.empty(0, dtype=np.intp), np.empty(0, dtype=np.intp), np.empty(0), np.empty(0))
array = np.asarray(result)
return array[:, 0].astype(np.intp), array[:, 1].astype(np.intp), array[:, 2], array[:, 3]
def _groups(sorted_keys: np.ndarray):
"""Yield (key, start, stop) for each run of equal values in a sorted array."""
if sorted_keys.size == 0:
return
boundaries = np.flatnonzero(sorted_keys[1:] != sorted_keys[:-1]) + 1
starts = np.concatenate(([0], boundaries))
stops = np.concatenate((boundaries, [sorted_keys.size]))
for start, stop in zip(starts.tolist(), stops.tolist(), strict=True):
yield int(sorted_keys[start]), start, stop
def _shell_components(vertices: np.ndarray, triangles: np.ndarray) -> list[np.ndarray]:
"""Edge-connected shells; weld identical source positions only for connectivity."""
triangles = np.asarray(triangles, dtype=np.intp).reshape(-1, 3)
if len(triangles) == 0:
return []
_, index = np.unique(vertices, axis=0, return_inverse=True)
faces = index.reshape(-1)[triangles]
edges = np.concatenate((faces[:, [0, 1]], faces[:, [1, 2]], faces[:, [2, 0]]), axis=0)
edges = np.sort(edges, axis=1)
owners = np.tile(np.arange(len(faces), dtype=np.intp), 3)
_, edge_ids = np.unique(edges, axis=0, return_inverse=True)
edge_ids = edge_ids.reshape(-1)
order = np.argsort(edge_ids, kind="stable")
sorted_ids = edge_ids[order]
sorted_owners = owners[order]
shared = sorted_ids[1:] == sorted_ids[:-1]
pairs = np.column_stack((sorted_owners[:-1][shared], sorted_owners[1:][shared]))
labels = np.arange(len(faces), dtype=np.intp)
while pairs.size:
lower = np.minimum(labels[pairs[:, 0]], labels[pairs[:, 1]])
updated = labels.copy()
np.minimum.at(updated, pairs[:, 0], lower)
np.minimum.at(updated, pairs[:, 1], lower)
updated = updated[updated] # Pointer jumping keeps the propagation logarithmic.
if np.array_equal(updated, labels):
break
labels = updated
# Each label is the smallest face index of its component, so this order is the
# order in which the components first appear.
return [triangles[np.flatnonzero(labels == label)] for label in np.unique(labels)]
def _shell_intervals(
vertices: np.ndarray,
triangles: np.ndarray,
east_axis: np.ndarray,
north_axis: np.ndarray,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""Intersect vertical grid rays with a closed triangle shell.
Return row, column, lower Z, upper Z for each interior interval. Multiple
intervals in one column retain openings and elevated structures. Triangle
edge duplicates are collapsed; tangencies have zero solid extent. Memory is
bounded by one building's intersected columns rather than a scene-wide cube.
"""
empty = (np.empty(0, dtype=np.intp), np.empty(0, dtype=np.intp), np.empty(0), np.empty(0))
vertices = np.asarray(vertices, dtype=np.float64)
triangles = np.asarray(triangles, dtype=np.intp).reshape(-1, 3)
east_axis = np.asarray(east_axis, dtype=np.float64)
north_axis = np.asarray(north_axis, dtype=np.float64)
width = len(east_axis)
if triangles.size == 0 or width == 0 or north_axis.size == 0:
return empty
corners = vertices[triangles]
a, b, c = corners[:, 0], corners[:, 1], corners[:, 2]
denominator = (b[:, 1] - c[:, 1]) * (a[:, 0] - c[:, 0]) + (c[:, 0] - b[:, 0]) * (
a[:, 1] - c[:, 1]
)
lower = corners.min(axis=1)
upper = corners.max(axis=1)
# Grid axes may run in either direction; candidate cells come from the sorted axes.
east_order = np.argsort(east_axis, kind="stable")
north_order = np.argsort(north_axis, kind="stable")
sorted_east = east_axis[east_order]
sorted_north = north_axis[north_order]
column_start = np.searchsorted(sorted_east, lower[:, 0], side="left")
column_stop = np.searchsorted(sorted_east, upper[:, 0], side="right")
row_start = np.searchsorted(sorted_north, lower[:, 1], side="left")
row_stop = np.searchsorted(sorted_north, upper[:, 1], side="right")
column_counts = np.maximum(column_stop - column_start, 0)
row_counts = np.maximum(row_stop - row_start, 0)
counts = column_counts * row_counts
counts[np.abs(denominator) < 1e-12] = 0 # Vertical walls do not cross a vertical ray.
candidates = np.flatnonzero(counts)
if candidates.size == 0:
return empty
keys: list[np.ndarray] = []
heights: list[np.ndarray] = []
for chunk in _chunks(candidates, counts[candidates]):
chunk_counts = counts[chunk]
total = int(chunk_counts.sum())
triangle_index = np.repeat(chunk, chunk_counts)
offsets = np.repeat(np.cumsum(chunk_counts) - chunk_counts, chunk_counts)
local = np.arange(total) - offsets
columns_per_row = np.repeat(column_counts[chunk], chunk_counts)
column_position = column_start[triangle_index] + local % columns_per_row
row_position = row_start[triangle_index] + local // columns_per_row
column_index = east_order[column_position]
row_index = north_order[row_position]
x = east_axis[column_index]
y = north_axis[row_index]
ta, tb, tc = a[triangle_index], b[triangle_index], c[triangle_index]
chunk_denominator = denominator[triangle_index]
u = (
(tb[:, 1] - tc[:, 1]) * (x - tc[:, 0]) + (tc[:, 0] - tb[:, 0]) * (y - tc[:, 1])
) / chunk_denominator
v = (
(tc[:, 1] - ta[:, 1]) * (x - tc[:, 0]) + (ta[:, 0] - tc[:, 0]) * (y - tc[:, 1])
) / chunk_denominator
inside = (u >= -1e-10) & (v >= -1e-10) & (u + v <= 1 + 1e-10)
z = u * ta[:, 2] + v * tb[:, 2] + (1 - u - v) * tc[:, 2]
keys.append(row_index[inside].astype(np.int64) * width + column_index[inside])
heights.append(z[inside])
key_array = np.concatenate(keys)
if key_array.size == 0:
return empty
z_array = np.concatenate(heights)
order = np.lexsort((z_array, key_array))
key_array = key_array[order]
z_array = z_array[order]
return _pair_levels(key_array, z_array, width)
def _chunks(candidates: np.ndarray, counts: np.ndarray):
"""Split triangles into runs whose candidate pair total stays within the memory bound."""
cumulative = np.cumsum(counts)
start = 0
while start < candidates.size:
limit = cumulative[start] - counts[start] + _PAIR_CHUNK
stop = int(np.searchsorted(cumulative, limit, side="right"))
stop = max(stop, start + 1)
yield candidates[start:stop]
start = stop
def _pair_levels(keys: np.ndarray, levels: np.ndarray, width: int):
"""Collapse near-duplicate crossings per column and pair them into intervals.
Columns without near-duplicate crossings are paired in one vectorized pass;
the few columns with crossings closer than the tolerance keep the sequential
rule that compares each value with the last retained level.
"""
same_column = keys[1:] == keys[:-1]
near = same_column & (levels[1:] - levels[:-1] <= _LEVEL_TOLERANCE)
irregular_keys = np.unique(keys[1:][near])
irregular = np.isin(keys, irregular_keys)
regular_keys = keys[~irregular]
regular_levels = levels[~irregular]
rows_out: list[np.ndarray] = []
lower_out: list[np.ndarray] = []
upper_out: list[np.ndarray] = []
if regular_keys.size:
boundaries = np.flatnonzero(regular_keys[1:] != regular_keys[:-1]) + 1
starts = np.concatenate(([0], boundaries))
stops = np.concatenate((boundaries, [regular_keys.size]))
sizes = stops - starts
odd = sizes % 2 == 1
if np.any(odd & (sizes > 1)):
raise ValueError("Native mesh has an unpaired vertical intersection; check topology")
paired = ~odd # Single crossings are boundary tangencies with zero occupied height.
pair_counts = sizes[paired] // 2
group_starts = starts[paired]
pair_offsets = np.repeat(np.cumsum(pair_counts) - pair_counts, pair_counts)
first = np.repeat(group_starts, pair_counts) + 2 * (
np.arange(pair_counts.sum()) - pair_offsets
)
rows_out.append(regular_keys[first])
lower_out.append(regular_levels[first])
upper_out.append(regular_levels[first + 1])
for key, start, stop in _groups(keys[irregular]):
merged: list[float] = []
for value in levels[irregular][start:stop].tolist():
if not merged or value - merged[-1] > _LEVEL_TOLERANCE:
merged.append(value)
if len(merged) % 2:
if len(merged) == 1:
continue
raise ValueError("Native mesh has an unpaired vertical intersection; check topology")
for lower, upper in zip(merged[::2], merged[1::2], strict=True):
if upper - lower > _LEVEL_TOLERANCE:
rows_out.append(np.asarray([key], dtype=np.int64))
lower_out.append(np.asarray([lower]))
upper_out.append(np.asarray([upper]))
if not rows_out:
return (np.empty(0, dtype=np.intp), np.empty(0, dtype=np.intp), np.empty(0), np.empty(0))
key_array = np.concatenate(rows_out)
lower_array = np.concatenate(lower_out)
upper_array = np.concatenate(upper_out)
order = np.lexsort((lower_array, key_array))
key_array = key_array[order]
return (
(key_array // width).astype(np.intp),
(key_array % width).astype(np.intp),
lower_array[order],
upper_array[order],
)
