"""CityJSON surface decoding without flattening roofs or filling courtyard holes."""
from __future__ import annotations
from collections.abc import Iterator
import numpy as np
import shapely
from shapely import constrained_delaunay_triangles
from shapely.geometry import Polygon
[docs]
def geometry_surfaces(geometry: dict) -> Iterator[tuple[list[list[int]], str]]:
"""Yield rings and matching semantic type at the CityJSON nesting depth."""
boundaries = geometry["boundaries"]
semantics = geometry.get("semantics", {})
values = semantics.get("values")
types = semantics.get("surfaces", [])
depth = {
"MultiSurface": 0,
"CompositeSurface": 0,
"Solid": 1,
"MultiSolid": 2,
"CompositeSolid": 2,
}.get(geometry["type"])
if depth is None:
raise ValueError(f"Unsupported CityJSON geometry: {geometry['type']}")
def visit(items, labels, remaining):
if labels is not None and len(labels) != len(items):
raise ValueError("CityJSON semantic nesting does not match boundaries")
for index, item in enumerate(items):
label = labels[index] if labels is not None else None
if remaining:
yield from visit(item, label, remaining - 1)
else:
if label is not None and (
not isinstance(label, int) or not 0 <= label < len(types)
):
raise ValueError("CityJSON semantic index is out of range")
yield item, types[label]["type"] if label is not None else "UnknownSurface"
yield from visit(boundaries, values, depth)
[docs]
def triangulate_surface(
vertices: np.ndarray,
rings: list[list[int]],
*,
planarity_tolerance_m: float = 0.05,
) -> list[tuple[int, int, int]]:
"""Constrained triangulation in the face plane, keeping original vertices/winding.
No synthetic roof height, vertex snapping, convex hull or fan triangulation
is used. Collinear zero-area faces produce no triangles. Non-planar,
self-intersecting or otherwise invalid faces are rejected for reporting.
"""
if not rings or any(len(ring) < 3 for ring in rings):
raise ValueError("CityJSON face has fewer than three vertices")
if any(index < 0 or index >= len(vertices) for ring in rings for index in ring):
raise ValueError("CityJSON vertex index is out of range")
points = np.asarray(vertices[rings[0]], dtype=np.float64)
centered = points - points[0]
# Newell normal also works for concave faces and vertical walls.
normal = _cross(centered, np.roll(centered, -1, axis=0)).sum(axis=0)
length = np.linalg.norm(normal)
if length < 1e-10:
# Source meshes can contain collapsed wall triangles. They enclose no
# area; dropping them must not also drop the rest of a valid building.
# A bow-tie ring also has a zero Newell normal, but is not collinear.
all_points = vertices[sorted({index for ring in rings for index in ring})]
singular = np.linalg.svd(all_points - all_points[0], compute_uv=False)
if len(singular) < 2 or singular[1] <= 1e-10:
return []
raise ValueError("CityJSON face is degenerate")
normal /= length
indices = sorted({index for ring in rings for index in ring})
if np.max(np.abs((vertices[indices] - points[0]) @ normal)) > planarity_tolerance_m:
raise ValueError("CityJSON face exceeds planarity tolerance")
if len(rings) == 1 and len(rings[0]) == 3:
return [tuple(rings[0])] # Native OBJ meshes are usually already triangulated.
keep = [axis for axis in range(3) if axis != int(np.argmax(np.abs(normal)))]
projected = vertices[:, keep]
polygon = Polygon(projected[rings[0]], [projected[ring] for ring in rings[1:]])
if not polygon.is_valid or polygon.area <= 0:
raise ValueError("CityJSON face has invalid rings")
lookup = {tuple(projected[index]): index for index in indices}
coordinates = shapely.get_coordinates(constrained_delaunay_triangles(polygon))
if coordinates.size == 0:
raise ValueError("CityJSON face produced no triangles")
corners = coordinates.reshape(-1, 4, 2)[:, :3, :].tolist()
triangles = np.asarray(
[[lookup[tuple(point)] for point in triangle] for triangle in corners], dtype=np.intp
)
a, b, c = (vertices[triangles[:, column]] for column in range(3))
facing = _cross(b - a, c - a)
# Summed in index order, as ``numpy.dot`` does for a vector of three.
against_normal = facing[:, 0] * normal[0] + facing[:, 1] * normal[1] + facing[:, 2] * normal[2]
triangles[against_normal < 0] = triangles[against_normal < 0][:, [0, 2, 1]]
return [tuple(triangle) for triangle in triangles.tolist()]
def _cross(first: np.ndarray, second: np.ndarray) -> np.ndarray:
"""Row-wise cross product with the component arithmetic of ``numpy.cross``."""
return np.column_stack(
(
first[:, 1] * second[:, 2] - first[:, 2] * second[:, 1],
first[:, 2] * second[:, 0] - first[:, 0] * second[:, 2],
first[:, 0] * second[:, 1] - first[:, 1] * second[:, 0],
)
)
