nodes/terrain_face.py

Part of Procedural Planets.

  1"""Single quad-sphere face: vectorised port of Lague's TerrainFace.
  2
  3`build_face_mesh()` produces a Mesh by:
  4
  51. Sampling an N×N grid of points on the unit cube along the given face axis.
  62. Normalising to the unit sphere.
  73. Asking ShapeGenerator for unscaled elevation (and updating its min/max).
  84. Displacing each vertex along its sphere normal by `radius * (1 + max(0, elev))`.
  95. Storing the *raw* elevation in UV.y (the shader / albedo lookup uses it later).
 106. Generating a triangle index list.
 11
 12Everything runs as whole-array numpy operations: one pass per face rather than
 13one per vertex. Run `harness.py` to time a six-face rebuild on your machine.
 14"""
 15
 16from __future__ import annotations
 17
 18import numpy as np
 19
 20from simvx.core.graphics.mesh import Mesh
 21
 22from .colour import BiomeColourSettings, biome_percent_array
 23from .shape import ShapeGenerator
 24
 25# Six face directions matching Lague's order (Vector3.up, down, left, right, forward, back).
 26# Y up, Y down, X-, X+, Z+, Z-.
 27FACE_DIRECTIONS: tuple[tuple[float, float, float], ...] = (
 28    (0.0, 1.0, 0.0),  # +Y
 29    (0.0, -1.0, 0.0),  # -Y
 30    (-1.0, 0.0, 0.0),  # -X
 31    (1.0, 0.0, 0.0),  # +X
 32    (0.0, 0.0, 1.0),  # +Z
 33    (0.0, 0.0, -1.0),  # -Z
 34)
 35
 36
 37def _axes_for(local_up: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
 38    """Reproduce Lague's axisA = (up.y, up.z, up.x); axisB = up × axisA."""
 39    axis_a = np.array([local_up[1], local_up[2], local_up[0]], dtype=np.float32)
 40    axis_b = np.cross(local_up, axis_a)
 41    return axis_a, axis_b
 42
 43
 44def build_face_indices(resolution: int) -> np.ndarray:
 45    """Index list for an N×N face grid: two triangles per cell, vectorised."""
 46    n = resolution
 47    # Grid of cells, each emits 2 triangles → 6 indices per cell.
 48    cell_count = (n - 1) * (n - 1)
 49    idx = np.empty(cell_count * 6, dtype=np.uint32)
 50
 51    # Top-left corner of every cell, flattened.
 52    yy, xx = np.meshgrid(np.arange(n - 1), np.arange(n - 1), indexing="ij")
 53    i = (xx + yy * n).ravel().astype(np.uint32)
 54
 55    # Lague's winding (Unity uses clockwise faces in left-handed):
 56    #   tri 0 = (i, i+n+1, i+n)
 57    #   tri 1 = (i, i+1,   i+n+1)
 58    idx[0::6] = i
 59    idx[1::6] = i + n + 1
 60    idx[2::6] = i + n
 61    idx[3::6] = i
 62    idx[4::6] = i + 1
 63    idx[5::6] = i + n + 1
 64    return idx
 65
 66
 67def build_face_mesh(
 68    local_up: tuple[float, float, float],
 69    resolution: int,
 70    shape: ShapeGenerator,
 71    biome_settings: BiomeColourSettings | None = None,
 72) -> Mesh:
 73    """Build a single face mesh and return it.
 74
 75    Updates ``shape.elevation_min`` / ``elevation_max`` as a side-effect.
 76    UV.x = biome % (0..1, sampled from biome_settings), drives the Y axis
 77    of the ramp texture. UV.y = raw elevation, drives the X axis. The two
 78    halves of the ramp (ocean / land) are selected by clamping elevation
 79    to [0, 1] in the consumer (see Planet.regenerate).
 80    """
 81    n = max(2, int(resolution))
 82    up = np.asarray(local_up, dtype=np.float32)
 83    axis_a, axis_b = _axes_for(up)
 84
 85    # Build the (N×N) grid of cube-surface points.
 86    coords = np.linspace(0.0, 1.0, n, dtype=np.float32) - 0.5  # percent - 0.5
 87    yy, xx = np.meshgrid(coords, coords, indexing="ij")
 88    # cube_pt = up + 2*xx*axis_a + 2*yy*axis_b  → shape (n, n, 3)
 89    cube_pts = (
 90        up[None, None, :]
 91        + (2.0 * xx)[..., None] * axis_a[None, None, :]
 92        + (2.0 * yy)[..., None] * axis_b[None, None, :]
 93    )
 94    cube_flat = cube_pts.reshape(-1, 3)
 95    # Normalise to unit sphere.
 96    lengths = np.linalg.norm(cube_flat, axis=1, keepdims=True)
 97    sphere_flat = cube_flat / np.maximum(lengths, 1e-8)
 98
 99    # Elevation pass: vectorised over all N² points.
100    raw_elev = shape.calculate_unscaled_elevation(sphere_flat)
101    scaled = shape.get_scaled_elevation(raw_elev)
102
103    # Displace along the sphere normal.
104    positions = sphere_flat * scaled[:, None]
105    normals = sphere_flat  # unit-sphere outward normals before displacement
106
107    # UVs: UV.x = biome %, UV.y = elevation %.
108    uvs = np.zeros((positions.shape[0], 2), dtype=np.float32)
109    if biome_settings is not None:
110        uvs[:, 0] = biome_percent_array(sphere_flat, biome_settings)
111    # Map raw elevation into [0, 1] using the running min/max: the renderer
112    # samples the X axis of the ramp from this. Clamp at the lower bound to
113    # 0.5 so anything ≤0 falls in the ocean half and >0 falls in the land half.
114    elev_min = shape.elevation_min if np.isfinite(shape.elevation_min) else 0.0
115    elev_max = shape.elevation_max if np.isfinite(shape.elevation_max) else 1.0
116    # Ocean (raw < 0) → x in [0, 0.5); land (raw >= 0) → x in [0.5, 1].
117    ocean_mask = raw_elev < 0.0
118    elev_uv = np.empty_like(raw_elev)
119    # Ocean half: map [elev_min, 0] -> [0.0, 0.5)
120    if ocean_mask.any():
121        ocean_lo = min(elev_min, 0.0)
122        ocean_pct = (raw_elev[ocean_mask] - ocean_lo) / max(0.0 - ocean_lo, 1e-3)
123        elev_uv[ocean_mask] = ocean_pct * 0.5
124    # Land half: map [0, elev_max] -> [0.5, 1.0]
125    land_mask = ~ocean_mask
126    if land_mask.any():
127        land_pct = raw_elev[land_mask] / max(elev_max, 1e-3)
128        elev_uv[land_mask] = 0.5 + np.clip(land_pct, 0.0, 1.0) * 0.5
129    uvs[:, 1] = elev_uv
130
131    indices = build_face_indices(n)
132
133    mesh = Mesh(
134        positions=positions.astype(np.float32),
135        indices=indices,
136        normals=normals.astype(np.float32),
137        texcoords=uvs,
138    )
139    return mesh