diff --git a/docs/api-reference.md b/docs/api-reference.md index a88dbc9..026cd02 100644 --- a/docs/api-reference.md +++ b/docs/api-reference.md @@ -350,6 +350,57 @@ Print and return memory breakdown of the current scene. --- +## LOD Utilities + +Level-of-detail helpers for terrain tiles and instanced geometry. + +```python +from rtxpy import compute_lod_level, compute_lod_distances, simplify_mesh, build_lod_chain +``` + +#### `compute_lod_level(distance, lod_distances)` + +Map a distance value to a discrete LOD level. + +| Parameter | Type | Description | +|-----------|------|-------------| +| `distance` | float | Distance from camera to object/tile center | +| `lod_distances` | list[float] | Ascending thresholds for LOD transitions | + +**Returns:** `int` — LOD level (0 = highest detail) + +#### `compute_lod_distances(tile_diagonal, factor=3.0, max_lod=3)` + +Generate distance thresholds from tile geometry. + +| Parameter | Type | Default | Description | +|-----------|------|---------|-------------| +| `tile_diagonal` | float | | Tile diagonal in world units | +| `factor` | float | `3.0` | Base multiplier for first threshold | +| `max_lod` | int | `3` | Number of LOD transitions | + +**Returns:** `list[float]` — distance thresholds + +#### `simplify_mesh(vertices, indices, ratio)` + +Simplify a triangle mesh via quadric decimation (requires `trimesh`). + +| Parameter | Type | Description | +|-----------|------|-------------| +| `vertices` | ndarray | Flat float32 vertex buffer `(N*3,)` | +| `indices` | ndarray | Flat int32 index buffer `(M*3,)` | +| `ratio` | float | Fraction of triangles to keep (0.0-1.0) | + +**Returns:** `(vertices, indices)` — simplified mesh buffers + +#### `build_lod_chain(vertices, indices, ratios=(1.0, 0.5, 0.25, 0.1))` + +Build a list of progressively simplified meshes. + +**Returns:** `list[(vertices, indices)]` — one pair per LOD level + +--- + ## RTX Class Low-level OptiX wrapper. Use this directly for custom ray tracing without the xarray accessor. diff --git a/docs/user-guide.md b/docs/user-guide.md index 4159e34..f763346 100644 --- a/docs/user-guide.md +++ b/docs/user-guide.md @@ -461,6 +461,7 @@ dem.rtx.explore(wind_data=wind) # Shift+W to toggle | **[** / **]** | Decrease / increase observer height | | **R** | Decrease terrain resolution (coarser) | | **Shift+R** | Increase terrain resolution (finer) | +| **Shift+A** | Toggle distance-based terrain LOD | | **Z** | Decrease vertical exaggeration | | **Shift+Z** | Increase vertical exaggeration | | **B** | Toggle mesh type (TIN / voxel) | @@ -537,6 +538,7 @@ write_stl('terrain.stl', verts, indices) ## Performance Tips +- **Enable terrain LOD**: Press `Shift+A` in the viewer to activate distance-based LOD. Nearby tiles render at full detail while distant tiles are automatically subsampled, cutting triangle count without visible quality loss - **Subsample large DEMs**: `dem[::2, ::2]` or `explore(subsample=4)` — 4x subsample is 16x less geometry - **Lower render_scale**: `explore(render_scale=0.25)` renders at quarter resolution for faster interaction - **Cache meshes**: Use `mesh_cache` parameter in `place_buildings()`, `place_roads()`, etc. to skip GeoJSON parsing on reload diff --git a/examples/terrain_lod_demo.py b/examples/terrain_lod_demo.py new file mode 100644 index 0000000..258c110 --- /dev/null +++ b/examples/terrain_lod_demo.py @@ -0,0 +1,100 @@ +"""Terrain LOD demo — distance-based level of detail for large terrains. + +Generates a synthetic 2048x2048 terrain and launches the interactive +viewer. Press Shift+A to toggle terrain LOD, which splits the terrain +into tiles and assigns each tile a resolution based on camera distance. + +Controls: + Shift+A Toggle terrain LOD on/off + R / Shift+R Manual resolution down / up (applies globally) + Z / Shift+Z Vertical exaggeration + WASD / arrows Move camera + +The LOD system is most useful with large terrains where full-resolution +rendering everywhere would exceed GPU memory or hurt frame rate. Nearby +tiles render at full detail while distant tiles are progressively +subsampled (2x, 4x, 8x). + +Usage: + python terrain_lod_demo.py + python terrain_lod_demo.py --size 4096 +""" + +import argparse + +import numpy as np +import xarray as xr + +import rtxpy # noqa: F401 — registers .rtx accessor + + +def make_terrain(size, seed=42): + """Generate a synthetic island terrain using multi-octave Perlin noise. + + Falls back to simple sine-based terrain if xarray-spatial is not + installed. + """ + try: + import cupy as cp + from xrspatial import generate_terrain + + template = xr.DataArray( + cp.zeros((size, size), dtype=cp.float32), dims=['y', 'x'], + ) + terrain = generate_terrain( + template, + x_range=(-5000, 5000), + y_range=(-5000, 5000), + seed=seed, + zfactor=3000, + full_extent=(-5000, -5000, 5000, 5000), + noise_mode='ridged', + warp_strength=0.3, + octaves=5, + ) + # Make it an island + sea_level = float(cp.nanmedian(terrain.data)) * 0.8 + terrain.data[:] = cp.maximum(terrain.data - sea_level, 0) + return terrain + except ImportError: + pass + + # Fallback: sin/cos-based terrain (no xarray-spatial needed) + rng = np.random.default_rng(seed) + y = np.linspace(0, 4 * np.pi, size, dtype=np.float32) + x = np.linspace(0, 4 * np.pi, size, dtype=np.float32) + yy, xx = np.meshgrid(y, x, indexing='ij') + z = (np.sin(xx) * np.cos(yy * 0.7) * 500 + + np.sin(xx * 2.3 + 1) * np.cos(yy * 1.7 + 0.5) * 200 + + rng.normal(0, 10, (size, size)).astype(np.float32)) + z = np.maximum(z, 0) + return xr.DataArray(z, dims=['y', 'x']) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser( + description="Terrain LOD demo — press Shift+A in viewer to toggle", + ) + parser.add_argument("--size", type=int, default=2048, + help="Grid size in pixels (default: 2048)") + parser.add_argument("--seed", type=int, default=42, + help="Random seed (default: 42)") + args = parser.parse_args() + + print(f"Generating {args.size}x{args.size} terrain...") + terrain = make_terrain(args.size, seed=args.seed) + + elev = terrain.data + if hasattr(elev, 'get'): + elev = elev.get() + print(f"Elevation range: {np.nanmin(elev):.0f} - {np.nanmax(elev):.0f} m") + print(f"\nPress Shift+A to toggle terrain LOD") + print(f"Press H for help overlay\n") + + terrain.rtx.explore( + width=1920, + height=1080, + render_scale=0.5, + subsample=1, + title="Terrain LOD Demo", + ) diff --git a/rtxpy/__init__.py b/rtxpy/__init__.py index ca3c137..a796878 100644 --- a/rtxpy/__init__.py +++ b/rtxpy/__init__.py @@ -23,6 +23,12 @@ load_meshes_from_zarr, chunks_for_pixel_window, ) +from .lod import ( + compute_lod_level, + compute_lod_distances, + simplify_mesh, + build_lod_chain, +) from .analysis import viewshed, hillshade, render, flyover, view from .engine import explore diff --git a/rtxpy/engine.py b/rtxpy/engine.py index fd7904b..642d3b4 100644 --- a/rtxpy/engine.py +++ b/rtxpy/engine.py @@ -739,9 +739,12 @@ def update(self, cam_x, cam_y, viewer): accessor._geometry_colors_dirty = True # Refresh viewer geometry tracking (same pattern as FIRMS toggle) + from .viewer.terrain_lod import is_terrain_lod_gid viewer._all_geometries = rtx.list_geometries() groups = set() for g in viewer._all_geometries: + if is_terrain_lod_gid(g): + continue parts = g.rsplit('_', 1) if len(parts) == 2 and parts[1].isdigit(): base = parts[0] @@ -1374,6 +1377,7 @@ class InteractiveViewer: - [/]: Decrease/increase observer height - R: Decrease terrain resolution (coarser, up to 8x subsample) - Shift+R: Increase terrain resolution (finer, down to 1x) + - Shift+A: Toggle distance-based terrain LOD - Z: Decrease vertical exaggeration - Shift+Z: Increase vertical exaggeration - B: Toggle mesh type (TIN / voxel) @@ -3235,6 +3239,22 @@ def _last_reload_time(self): def _last_reload_time(self, value): self.terrain._last_reload_time = value + @property + def lod_enabled(self): + return self.terrain.lod_enabled + + @lod_enabled.setter + def lod_enabled(self, value): + self.terrain.lod_enabled = value + + @property + def _terrain_lod_manager(self): + return self.terrain._terrain_lod_manager + + @_terrain_lod_manager.setter + def _terrain_lod_manager(self, value): + self.terrain._terrain_lod_manager = value + # ------------------------------------------------------------------ # Delegation properties — GeometryLayerManager # ------------------------------------------------------------------ @@ -3615,6 +3635,49 @@ def _apply_cyberpunk_minimap(self, rgba, water): result[:, :, :3] = np.clip(result[:, :, :3], 0, 1) return result + def _enable_terrain_lod(self): + """Switch terrain rendering from single-GAS to per-tile LOD. + + Removes the single ``'terrain'`` (and ``'terrain_skirt'``) + geometry and creates a :class:`TerrainLODManager` that renders + each tile at a distance-appropriate resolution. + """ + from .viewer.terrain_lod import TerrainLODManager + + # Get full-res terrain as numpy + base = self._base_raster + terrain_data = base.data + if hasattr(terrain_data, 'get'): + terrain_np = terrain_data.get() + else: + terrain_np = np.asarray(terrain_data) + + # Remove the single terrain geometry + if self.rtx.has_geometry('terrain'): + self.rtx.remove_geometry('terrain') + if self.rtx.has_geometry('terrain_skirt'): + self.rtx.remove_geometry('terrain_skirt') + + # Choose tile size: aim for ~8-16 tiles across largest dimension + H, W = terrain_np.shape + tile_size = max(32, min(256, max(H, W) // 8)) + + mgr = TerrainLODManager( + terrain_np, + tile_size=tile_size, + pixel_spacing_x=self._base_pixel_spacing_x, + pixel_spacing_y=self._base_pixel_spacing_y, + max_lod=3, + base_subsample=self.subsample_factor, + ) + self._terrain_lod_manager = mgr + self.lod_enabled = True + + # Force initial tile build + mgr.update(self.position, self.rtx, + ve=self.vertical_exaggeration, force=True) + self._update_frame() + def _rebuild_at_resolution(self, factor): """Rebuild terrain mesh at a different subsample factor. @@ -3630,6 +3693,17 @@ def _rebuild_at_resolution(self, factor): from . import mesh as mesh_mod self.subsample_factor = factor + + # If LOD is active, update the LOD manager's base subsample + # and force a full tile rebuild instead of the single-GAS path. + if self.lod_enabled and self._terrain_lod_manager is not None: + self._terrain_lod_manager.set_base_subsample(factor) + self._terrain_lod_manager.update( + self.position, self.rtx, + ve=self.vertical_exaggeration, force=True) + # Still need to update raster/spacing for overlays and re-snapping + # but skip the single-terrain GAS rebuild below. + base = self._base_raster # 1. Subsample the raster @@ -3658,7 +3732,16 @@ def _rebuild_at_resolution(self, factor): ve = self.vertical_exaggeration cache_key = (factor, self.mesh_type) - if self.mesh_type == 'heightfield': + # When LOD is active, tiles are built by the LOD manager. + # We still need terrain_np for elevation stats (computed below). + _lod_active = (self.lod_enabled and self._terrain_lod_manager is not None) + if _lod_active: + terrain_data = sub.data + if hasattr(terrain_data, 'get'): + terrain_np = terrain_data.get() + else: + terrain_np = np.asarray(terrain_data) + elif self.mesh_type == 'heightfield': # Heightfield path: no triangle mesh needed if cache_key in self._terrain_mesh_cache: _, _, terrain_np = self._terrain_mesh_cache[cache_key] @@ -3801,8 +3884,11 @@ def _rebuild_at_resolution(self, factor): if has_cupy: gpu_terrain = cp.asarray(terrain_np) self._gpu_terrain = gpu_terrain + from .viewer.terrain_lod import is_terrain_lod_gid for geom_id in self.rtx.list_geometries(): - if geom_id == 'terrain': + if geom_id == 'terrain' or geom_id == 'terrain_skirt': + continue + if is_terrain_lod_gid(geom_id): continue # Baked meshes — re-snap Z to new terrain surface + VE if hasattr(self, '_baked_meshes') and geom_id in self._baked_meshes: @@ -3923,10 +4009,19 @@ def _rebuild_vertical_exaggeration(self, ve): self.vertical_exaggeration = ve H, W = self.terrain_shape - # Use cached mesh if available, otherwise build and cache - cache_key = (self.subsample_factor, self.mesh_type) - - if self.mesh_type == 'heightfield': + # If LOD is active, force a full tile rebuild with the new VE + if self.lod_enabled and self._terrain_lod_manager is not None: + self._terrain_lod_manager._tile_cache.clear() + self._terrain_lod_manager._tile_lods.clear() + self._terrain_lod_manager.update( + self.position, self.rtx, ve=ve, force=True) + # Still need terrain_np for elevation stats below + terrain_data = self.raster.data + if hasattr(terrain_data, 'get'): + terrain_np = terrain_data.get() + else: + terrain_np = np.asarray(terrain_data) + elif self.mesh_type == 'heightfield': # Heightfield path: rebuild GAS with new VE if cache_key in self._terrain_mesh_cache: _, _, terrain_np = self._terrain_mesh_cache[cache_key] @@ -4026,8 +4121,11 @@ def _rebuild_vertical_exaggeration(self, ve): if has_cupy: gpu_terrain = cp.asarray(terrain_np) self._gpu_terrain = gpu_terrain + from .viewer.terrain_lod import is_terrain_lod_gid for geom_id in self.rtx.list_geometries(): - if geom_id == 'terrain': + if geom_id == 'terrain' or geom_id == 'terrain_skirt': + continue + if is_terrain_lod_gid(geom_id): continue # Baked meshes (merged buildings/curves) — re-snap Z to terrain + VE if hasattr(self, '_baked_meshes') and geom_id in self._baked_meshes: @@ -4464,9 +4562,12 @@ def _toggle_firms(self): # Refresh geometry layer tracking if self.rtx is not None: + from .viewer.terrain_lod import is_terrain_lod_gid self._all_geometries = self.rtx.list_geometries() groups = set() for g in self._all_geometries: + if is_terrain_lod_gid(g): + continue parts = g.rsplit('_', 1) if len(parts) == 2 and parts[1].isdigit(): base = parts[0] @@ -6356,12 +6457,27 @@ def _action_toggle_wind(self): self._toggle_wind() def _action_toggle_terrain_vis(self): + from .viewer.terrain_lod import is_terrain_lod_gid entry = self.rtx._geom_state.gas_entries.get('terrain') if entry is not None: vis = not entry.visible self.rtx.set_geometry_visible('terrain', vis) print(f"Terrain {'shown' if vis else 'hidden'}") self._needs_render = True + elif self.lod_enabled: + # Determine current visibility from any LOD tile + vis = True + for gid in self.rtx.list_geometries(): + if is_terrain_lod_gid(gid): + e = self.rtx._geom_state.gas_entries.get(gid) + if e is not None: + vis = not e.visible + break + for gid in self.rtx.list_geometries(): + if is_terrain_lod_gid(gid): + self.rtx.set_geometry_visible(gid, vis) + print(f"Terrain {'shown' if vis else 'hidden'}") + self._needs_render = True def _action_toggle_gtfs_rt(self): self._toggle_gtfs_rt() @@ -6492,6 +6608,25 @@ def _action_resolution_finer(self): if new_factor != self.subsample_factor: self._rebuild_at_resolution(new_factor) + def _action_toggle_terrain_lod(self): + """Toggle distance-based terrain LOD on/off.""" + if self.rtx is None: + return + + if self.lod_enabled: + # Disable LOD — remove tile geometries, restore single terrain + if self._terrain_lod_manager is not None: + self._terrain_lod_manager.remove_all(self.rtx) + self._terrain_lod_manager = None + self.lod_enabled = False + # Rebuild the single terrain geometry + self._rebuild_at_resolution(self.subsample_factor) + print("Terrain LOD: OFF") + else: + # Enable LOD — replace single terrain with tiled LOD + self._enable_terrain_lod() + print(f"Terrain LOD: ON ({self._terrain_lod_manager.get_stats()})") + def _action_ve_down(self): new_ve = max(0.1, round(self.vertical_exaggeration - 0.1, 1)) if new_ve != self.vertical_exaggeration: @@ -6984,6 +7119,12 @@ def _tick(self): if self._chunk_manager.update(self.position[0], self.position[1], self): self._geometry_colors_builder = self._accessor._build_geometry_colors_gpu self._render_needed = True + # Terrain LOD: update tile resolutions based on camera distance + if self.lod_enabled and self._terrain_lod_manager is not None: + if self._terrain_lod_manager.update( + self.position, self.rtx, + ve=self.vertical_exaggeration): + self._render_needed = True # AO/DOF: keep accumulating samples when camera is stationary if ((self.ao_enabled or self.dof_enabled) and not self._held_keys and not self._mouse_dragging @@ -8831,6 +8972,7 @@ def _render_help_text(self): ("Y", "Cycle color stretch"), (", / .", "Overlay alpha"), ("R / Shift+R", "Resolution down / up"), + ("Shift+A", "Toggle terrain LOD"), ("Z / Shift+Z", "Vert. exag. down / up"), ("B", "Toggle TIN / Voxel"), ("T", "Toggle shadows"), @@ -9480,6 +9622,7 @@ def explore(raster, width: int = 800, height: int = 600, - [/]: Decrease/increase observer height - R: Decrease terrain resolution (coarser, up to 8x subsample) - Shift+R: Increase terrain resolution (finer, down to 1x) + - Shift+A: Toggle distance-based terrain LOD - Z: Decrease vertical exaggeration - Shift+Z: Increase vertical exaggeration - B: Toggle mesh type (TIN / voxel) diff --git a/rtxpy/lod.py b/rtxpy/lod.py new file mode 100644 index 0000000..c25d050 --- /dev/null +++ b/rtxpy/lod.py @@ -0,0 +1,131 @@ +"""Level-of-detail utilities for terrain and instanced geometry. + +Provides LOD level computation, mesh simplification via quadric +decimation, and LOD chain generation for multi-resolution rendering. +""" + +import numpy as np + + +def compute_lod_level(distance, lod_distances): + """Return the LOD level for a given distance. + + Parameters + ---------- + distance : float + Distance from camera to object or tile center. + lod_distances : list of float + Ascending distance thresholds defining LOD transitions. + For *N* thresholds, LOD levels range from 0 to *N*. + Level 0 is the highest detail (closest to camera). + + Returns + ------- + int + LOD level (0 = highest detail). + + Examples + -------- + >>> compute_lod_level(100, [500, 1000, 2000]) + 0 + >>> compute_lod_level(750, [500, 1000, 2000]) + 1 + >>> compute_lod_level(3000, [500, 1000, 2000]) + 3 + """ + for i, threshold in enumerate(lod_distances): + if distance < threshold: + return i + return len(lod_distances) + + +def compute_lod_distances(tile_diagonal, factor=3.0, max_lod=3): + """Compute LOD distance thresholds from tile geometry. + + Each LOD transition occurs at ``tile_diagonal * factor * 2^level``. + + Parameters + ---------- + tile_diagonal : float + Diagonal size of a terrain tile in world units. + factor : float + Base multiplier for the first threshold. + max_lod : int + Number of LOD transitions (max LOD level = max_lod). + + Returns + ------- + list of float + Distance thresholds for LOD 0 → 1, 1 → 2, etc. + """ + return [tile_diagonal * factor * (2 ** i) for i in range(max_lod)] + + +def simplify_mesh(vertices, indices, ratio): + """Simplify a triangle mesh using quadric decimation. + + Requires ``trimesh`` with ``simplify_quadric_decimation`` support. + Falls back to returning the original mesh if trimesh is unavailable + or simplification fails. + + Parameters + ---------- + vertices : np.ndarray + Flat float32 vertex buffer, shape ``(N*3,)``. + indices : np.ndarray + Flat int32 index buffer, shape ``(M*3,)``. + ratio : float + Fraction of triangles to keep (0.0 to 1.0). + + Returns + ------- + simplified_vertices : np.ndarray + Flat float32 vertex buffer. + simplified_indices : np.ndarray + Flat int32 index buffer. + """ + if ratio >= 1.0: + return vertices, indices + + try: + import trimesh + except ImportError: + return vertices, indices + + verts_2d = np.asarray(vertices, dtype=np.float32).reshape(-1, 3) + faces_2d = np.asarray(indices, dtype=np.int32).reshape(-1, 3) + num_target = max(4, int(len(faces_2d) * ratio)) + + try: + mesh = trimesh.Trimesh(vertices=verts_2d, faces=faces_2d, + process=False) + simplified = mesh.simplify_quadric_decimation(num_target) + return (simplified.vertices.astype(np.float32).flatten(), + simplified.faces.astype(np.int32).flatten()) + except Exception: + return vertices, indices + + +def build_lod_chain(vertices, indices, ratios=(1.0, 0.5, 0.25, 0.1)): + """Build a chain of progressively simplified meshes. + + Parameters + ---------- + vertices : np.ndarray + Full-detail flat float32 vertex buffer. + indices : np.ndarray + Full-detail flat int32 index buffer. + ratios : tuple of float + Triangle retention ratios per LOD level. The first entry + should be 1.0 (original mesh). + + Returns + ------- + list of (np.ndarray, np.ndarray) + ``(vertices, indices)`` tuples, one per LOD level. + """ + chain = [(vertices.copy(), indices.copy())] + for ratio in ratios[1:]: + v, i = simplify_mesh(vertices, indices, ratio) + chain.append((v, i)) + return chain diff --git a/rtxpy/tests/test_lod.py b/rtxpy/tests/test_lod.py new file mode 100644 index 0000000..8147a82 --- /dev/null +++ b/rtxpy/tests/test_lod.py @@ -0,0 +1,333 @@ +"""Tests for level-of-detail utilities and terrain LOD manager.""" + +import numpy as np +import pytest + +from rtxpy.lod import ( + compute_lod_level, + compute_lod_distances, + simplify_mesh, + build_lod_chain, +) +from rtxpy.viewer.terrain_lod import ( + TerrainLODManager, + is_terrain_lod_gid, + _add_tile_skirt, + _tile_gid, +) + + +# --------------------------------------------------------------------------- +# compute_lod_level +# --------------------------------------------------------------------------- + +class TestComputeLodLevel: + """Tests for distance-to-LOD mapping.""" + + def test_below_first_threshold(self): + assert compute_lod_level(10, [100, 200, 400]) == 0 + + def test_between_thresholds(self): + assert compute_lod_level(150, [100, 200, 400]) == 1 + + def test_above_all_thresholds(self): + assert compute_lod_level(500, [100, 200, 400]) == 3 + + def test_exact_boundary(self): + # At the exact boundary, distance < threshold is false → next level + assert compute_lod_level(100, [100, 200, 400]) == 1 + + def test_single_threshold(self): + assert compute_lod_level(50, [100]) == 0 + assert compute_lod_level(150, [100]) == 1 + + def test_empty_thresholds(self): + assert compute_lod_level(50, []) == 0 + + def test_zero_distance(self): + assert compute_lod_level(0, [100, 200]) == 0 + + def test_negative_distance(self): + assert compute_lod_level(-10, [100, 200]) == 0 + + +# --------------------------------------------------------------------------- +# compute_lod_distances +# --------------------------------------------------------------------------- + +class TestComputeLodDistances: + """Tests for LOD distance threshold generation.""" + + def test_basic(self): + dists = compute_lod_distances(100.0, factor=2.0, max_lod=3) + assert len(dists) == 3 + assert dists[0] == pytest.approx(200.0) + assert dists[1] == pytest.approx(400.0) + assert dists[2] == pytest.approx(800.0) + + def test_single_level(self): + dists = compute_lod_distances(50.0, factor=3.0, max_lod=1) + assert len(dists) == 1 + assert dists[0] == pytest.approx(150.0) + + def test_zero_max_lod(self): + dists = compute_lod_distances(100.0, factor=2.0, max_lod=0) + assert dists == [] + + +# --------------------------------------------------------------------------- +# simplify_mesh +# --------------------------------------------------------------------------- + +class TestSimplifyMesh: + """Tests for mesh simplification.""" + + def _make_grid_mesh(self, H=10, W=10): + """Create a simple grid mesh for testing.""" + n_verts = H * W + n_tris = (H - 1) * (W - 1) * 2 + verts = np.zeros(n_verts * 3, dtype=np.float32) + indices = np.zeros(n_tris * 3, dtype=np.int32) + + for h in range(H): + for w in range(W): + idx = (h * W + w) * 3 + verts[idx] = float(w) + verts[idx + 1] = float(h) + verts[idx + 2] = float(h + w) * 0.1 + + tri = 0 + for h in range(H - 1): + for w in range(W - 1): + v00 = h * W + w + v01 = h * W + w + 1 + v10 = (h + 1) * W + w + v11 = (h + 1) * W + w + 1 + indices[tri * 3] = v00 + indices[tri * 3 + 1] = v11 + indices[tri * 3 + 2] = v10 + tri += 1 + indices[tri * 3] = v00 + indices[tri * 3 + 1] = v01 + indices[tri * 3 + 2] = v11 + tri += 1 + + return verts, indices + + def test_ratio_one_returns_original(self): + verts, indices = self._make_grid_mesh() + sv, si = simplify_mesh(verts, indices, 1.0) + np.testing.assert_array_equal(sv, verts) + np.testing.assert_array_equal(si, indices) + + def test_simplification_reduces_triangles(self): + verts, indices = self._make_grid_mesh(20, 20) + orig_tris = len(indices) // 3 + sv, si = simplify_mesh(verts, indices, 0.5) + new_tris = len(si) // 3 + # Simplified mesh should have fewer or equal triangles + assert new_tris <= orig_tris + + +# --------------------------------------------------------------------------- +# build_lod_chain +# --------------------------------------------------------------------------- + +class TestBuildLodChain: + def test_chain_length(self): + verts = np.zeros(30, dtype=np.float32) + indices = np.zeros(12, dtype=np.int32) + chain = build_lod_chain(verts, indices, ratios=(1.0, 0.5, 0.25)) + assert len(chain) == 3 + + def test_first_level_is_copy(self): + verts = np.arange(12, dtype=np.float32) + indices = np.arange(6, dtype=np.int32) + chain = build_lod_chain(verts, indices, ratios=(1.0,)) + np.testing.assert_array_equal(chain[0][0], verts) + np.testing.assert_array_equal(chain[0][1], indices) + # Must be a copy, not the same object + assert chain[0][0] is not verts + + +# --------------------------------------------------------------------------- +# TerrainLODManager +# --------------------------------------------------------------------------- + +class _FakeRTX: + """Minimal stub for RTX in tests — tracks add/remove calls.""" + + def __init__(self): + self.geometries = {} + + def add_geometry(self, gid, verts, indices, **kw): + self.geometries[gid] = (verts.copy(), indices.copy()) + return 0 + + def remove_geometry(self, gid): + self.geometries.pop(gid, None) + return 0 + + def has_geometry(self, gid): + return gid in self.geometries + + def list_geometries(self): + return list(self.geometries.keys()) + + +class TestTerrainLODManager: + """Tests for the terrain LOD tile manager.""" + + def _make_terrain(self, H=256, W=256): + """Create a synthetic terrain: gentle gradient.""" + y = np.linspace(0, 100, H, dtype=np.float32) + x = np.linspace(0, 100, W, dtype=np.float32) + return y[:, None] + x[None, :] + + def test_tile_count(self): + terrain = self._make_terrain(256, 256) + mgr = TerrainLODManager(terrain, tile_size=64) + assert mgr.n_tiles == 16 # 4x4 + + def test_tile_count_non_divisible(self): + terrain = self._make_terrain(300, 200) + mgr = TerrainLODManager(terrain, tile_size=128) + # 300/128 = ceil(2.34) = 3 rows, 200/128 = ceil(1.56) = 2 cols + assert mgr.n_tiles == 6 + + def test_initial_update_creates_tiles(self): + terrain = self._make_terrain(128, 128) + rtx = _FakeRTX() + mgr = TerrainLODManager(terrain, tile_size=64, + pixel_spacing_x=1.0, pixel_spacing_y=1.0) + changed = mgr.update(np.array([64, 64, 0]), rtx, force=True) + assert changed + # Should have created 4 tiles (2x2 grid) + assert mgr.n_tiles == 4 + for tr in range(2): + for tc in range(2): + gid = _tile_gid(tr, tc) + assert gid in rtx.geometries, f"Missing tile {gid}" + + def test_no_change_without_movement(self): + terrain = self._make_terrain(128, 128) + rtx = _FakeRTX() + mgr = TerrainLODManager(terrain, tile_size=64) + mgr.update(np.array([64, 64, 0]), rtx, force=True) + changed = mgr.update(np.array([64, 64, 0]), rtx) + assert not changed + + def test_lod_varies_with_distance(self): + terrain = self._make_terrain(512, 512) + rtx = _FakeRTX() + mgr = TerrainLODManager( + terrain, tile_size=128, + pixel_spacing_x=1.0, pixel_spacing_y=1.0, + max_lod=3, lod_distance_factor=1.0, + ) + # Camera at corner — near tiles get LOD 0, far tiles get higher LOD + mgr.update(np.array([0, 0, 0]), rtx, force=True) + lods = mgr.tile_lods + # Tile (0,0) is nearest, should be LOD 0 + assert lods[(0, 0)] == 0 + # Tile (3,3) is farthest, should have higher LOD + assert lods[(3, 3)] > lods[(0, 0)] + + def test_remove_all_clears_scene(self): + terrain = self._make_terrain(128, 128) + rtx = _FakeRTX() + mgr = TerrainLODManager(terrain, tile_size=64) + mgr.update(np.array([64, 64, 0]), rtx, force=True) + assert len(rtx.geometries) > 0 + mgr.remove_all(rtx) + # All terrain LOD tiles should be removed + for gid in rtx.geometries: + assert not is_terrain_lod_gid(gid) + + def test_ve_applied_to_z(self): + terrain = self._make_terrain(64, 64) + rtx = _FakeRTX() + mgr = TerrainLODManager(terrain, tile_size=64) + mgr.update(np.array([32, 32, 0]), rtx, ve=2.0, force=True) + gid = _tile_gid(0, 0) + verts = rtx.geometries[gid][0] + z_vals = verts[2::3] + # With VE=2, z values should be doubled relative to terrain + assert float(np.max(z_vals)) > float(np.max(terrain)) + + def test_base_subsample_changes_invalidate_cache(self): + terrain = self._make_terrain(128, 128) + rtx = _FakeRTX() + mgr = TerrainLODManager(terrain, tile_size=64, base_subsample=1) + mgr.update(np.array([64, 64, 0]), rtx, force=True) + gid = _tile_gid(0, 0) + verts1 = rtx.geometries[gid][0].copy() + + mgr.set_base_subsample(2) + mgr.update(np.array([64, 64, 0]), rtx, force=True) + verts2 = rtx.geometries[gid][0] + # Meshes should differ (fewer vertices at 2x subsample) + assert len(verts2) != len(verts1) + + def test_get_stats(self): + terrain = self._make_terrain(128, 128) + mgr = TerrainLODManager(terrain, tile_size=64) + # Before update + stats = mgr.get_stats() + assert "no tiles" in stats + # After update + rtx = _FakeRTX() + mgr.update(np.array([64, 64, 0]), rtx, force=True) + stats = mgr.get_stats() + assert "LOD tiles:" in stats + + +# --------------------------------------------------------------------------- +# Tile helpers +# --------------------------------------------------------------------------- + +class TestTileHelpers: + def test_tile_gid_format(self): + assert _tile_gid(0, 0) == "terrain_lod_r0_c0" + assert _tile_gid(3, 7) == "terrain_lod_r3_c7" + + def test_is_terrain_lod_gid(self): + assert is_terrain_lod_gid("terrain_lod_r0_c0") + assert is_terrain_lod_gid("terrain_lod_r12_c34") + assert not is_terrain_lod_gid("terrain") + assert not is_terrain_lod_gid("terrain_skirt") + assert not is_terrain_lod_gid("buildings_0") + + +class TestAddTileSkirt: + def test_adds_skirt_vertices(self): + """Skirt should add perimeter + wall vertices.""" + H, W = 4, 4 + n_verts = H * W + n_tris = (H - 1) * (W - 1) * 2 + + verts = np.zeros(n_verts * 3, dtype=np.float32) + indices = np.zeros(n_tris * 3, dtype=np.int32) + for h in range(H): + for w in range(W): + idx = (h * W + w) * 3 + verts[idx] = float(w) + verts[idx + 1] = float(h) + verts[idx + 2] = float(h + w) + + new_v, new_i = _add_tile_skirt(verts, indices, H, W) + # Perimeter of 4x4 grid: 4+3+3+2 = 12 vertices added + n_perim = 2 * (H + W) - 4 + assert len(new_v) == (n_verts + n_perim) * 3 + assert len(new_i) > len(indices) + + def test_skirt_z_below_min(self): + H, W = 3, 3 + verts = np.zeros(9 * 3, dtype=np.float32) + for i in range(9): + verts[i * 3 + 2] = 10.0 # all z = 10 + indices = np.zeros(8 * 3, dtype=np.int32) + + new_v, _ = _add_tile_skirt(verts, indices, H, W) + skirt_z = new_v[9 * 3 + 2::3] # z of skirt vertices + assert np.all(skirt_z < 10.0) diff --git a/rtxpy/viewer/keybindings.py b/rtxpy/viewer/keybindings.py index 4d0aefa..82d1c93 100644 --- a/rtxpy/viewer/keybindings.py +++ b/rtxpy/viewer/keybindings.py @@ -40,6 +40,7 @@ 'T': '_action_cycle_time', # Time-of-day 'Y': '_action_toggle_hydro', # Hydro flow particles 'N': '_action_toggle_clouds', # Cloud layer + 'A': '_action_toggle_terrain_lod', # Distance-based terrain LOD } # Lowercase key bindings — checked after shift bindings diff --git a/rtxpy/viewer/terrain.py b/rtxpy/viewer/terrain.py index c57e1b3..9e8064d 100644 --- a/rtxpy/viewer/terrain.py +++ b/rtxpy/viewer/terrain.py @@ -24,6 +24,8 @@ class TerrainState: '_coord_step_x', '_coord_step_y', '_reload_cooldown', '_last_reload_time', '_terrain_reload_future', '_terrain_reload_pool', + # LOD state + 'lod_enabled', '_terrain_lod_manager', ) def __init__(self, raster, pixel_spacing_x=1.0, pixel_spacing_y=1.0, @@ -63,3 +65,7 @@ def __init__(self, raster, pixel_spacing_x=1.0, pixel_spacing_y=1.0, self._last_reload_time = 0.0 self._terrain_reload_future = None self._terrain_reload_pool = None + + # LOD + self.lod_enabled = False + self._terrain_lod_manager = None diff --git a/rtxpy/viewer/terrain_lod.py b/rtxpy/viewer/terrain_lod.py new file mode 100644 index 0000000..00b852d --- /dev/null +++ b/rtxpy/viewer/terrain_lod.py @@ -0,0 +1,321 @@ +"""Distance-based terrain LOD manager. + +Divides a terrain raster into a grid of tiles. Each tile is assigned +an LOD level based on its distance to the camera, controlling the mesh +resolution (LOD 0 = full detail, LOD *k* = 2^k subsampling). Tiles +are built as individual GAS entries in the OptiX IAS so the raytracer +traverses only the detail actually needed. + +Tile edges get a short vertical skirt to hide T-junction cracks where +adjacent tiles have different LOD levels. +""" + +import numpy as np + +from ..lod import compute_lod_level, compute_lod_distances + + +class TerrainLODManager: + """Manages per-tile terrain LOD in the interactive viewer. + + Parameters + ---------- + terrain_np : np.ndarray + Full-resolution elevation array, shape ``(H, W)``. + tile_size : int + Tile edge length in full-resolution pixels (default 128). + pixel_spacing_x, pixel_spacing_y : float + World-space size of one raster pixel. + max_lod : int + Maximum LOD level. Level *k* subsamples the tile by ``2^k``. + lod_distance_factor : float + Controls how far each LOD band extends. First transition + at ``tile_diagonal * factor``, doubling per level. + base_subsample : int + Global base subsample factor (from R/Shift+R). LOD 0 uses + this value; LOD *k* uses ``base_subsample * 2^k``. + """ + + __slots__ = ( + '_terrain_np', '_tile_size', '_psx', '_psy', + '_max_lod', '_lod_distances', '_lod_distance_factor', + '_H', '_W', '_n_tile_rows', '_n_tile_cols', + '_tile_centers', '_tile_lods', '_tile_cache', + '_active_tiles', '_last_update_pos', '_update_threshold', + '_base_subsample', + ) + + def __init__(self, terrain_np, tile_size=128, + pixel_spacing_x=1.0, pixel_spacing_y=1.0, + max_lod=3, lod_distance_factor=3.0, + base_subsample=1): + self._terrain_np = terrain_np + self._tile_size = tile_size + self._psx = pixel_spacing_x + self._psy = pixel_spacing_y + self._max_lod = max_lod + self._lod_distance_factor = lod_distance_factor + self._base_subsample = base_subsample + + H, W = terrain_np.shape + self._H = H + self._W = W + self._n_tile_rows = (H + tile_size - 1) // tile_size + self._n_tile_cols = (W + tile_size - 1) // tile_size + + # Pre-compute tile centres in world coordinates + self._tile_centers = {} + for tr in range(self._n_tile_rows): + for tc in range(self._n_tile_cols): + r0 = tr * tile_size + c0 = tc * tile_size + r1 = min(r0 + tile_size, H) + c1 = min(c0 + tile_size, W) + cx = (c0 + c1) * 0.5 * pixel_spacing_x + cy = (r0 + r1) * 0.5 * pixel_spacing_y + self._tile_centers[(tr, tc)] = (cx, cy) + + # LOD distance thresholds + tile_diag = np.sqrt( + (tile_size * pixel_spacing_x) ** 2 + + (tile_size * pixel_spacing_y) ** 2 + ) + self._lod_distances = compute_lod_distances( + tile_diag, factor=lod_distance_factor, max_lod=max_lod) + + # Per-tile state + self._tile_lods = {} # (tr, tc) -> current LOD level + self._tile_cache = {} # (tr, tc, lod, base_sub) -> (verts, indices) + self._active_tiles = set() # GAS IDs currently in the scene + + # Movement threshold before re-evaluating LOD + self._last_update_pos = None + self._update_threshold = tile_diag * 0.25 + + # ------------------------------------------------------------------ + # Public API + # ------------------------------------------------------------------ + + @property + def n_tiles(self): + """Total number of terrain tiles.""" + return self._n_tile_rows * self._n_tile_cols + + @property + def tile_lods(self): + """Current LOD assignment per tile: ``{(row, col): level}``.""" + return dict(self._tile_lods) + + def set_base_subsample(self, factor): + """Update the global base subsample and invalidate cache.""" + if factor != self._base_subsample: + self._base_subsample = factor + self._tile_cache.clear() + self._tile_lods.clear() + + def set_terrain(self, terrain_np): + """Replace the terrain data (e.g. after dynamic reload).""" + self._terrain_np = terrain_np + H, W = terrain_np.shape + self._H = H + self._W = W + self._tile_cache.clear() + self._tile_lods.clear() + + def update(self, camera_pos, rtx, ve=1.0, force=False): + """Re-evaluate LOD per tile and rebuild changed tiles. + + Parameters + ---------- + camera_pos : array-like + Camera position ``[x, y, z]`` in world coordinates. + rtx : RTX + Scene handle for adding/removing tile GAS. + ve : float + Current vertical exaggeration. + force : bool + If True, rebuild all tiles regardless of movement threshold. + + Returns + ------- + bool + True if any tile GAS was added or updated. + """ + cam_x, cam_y = float(camera_pos[0]), float(camera_pos[1]) + + # Skip if camera hasn't moved enough + if not force and self._last_update_pos is not None: + dx = cam_x - self._last_update_pos[0] + dy = cam_y - self._last_update_pos[1] + if dx * dx + dy * dy < self._update_threshold ** 2: + return False + + self._last_update_pos = (cam_x, cam_y) + + changed = False + new_tile_ids = set() + + for tr in range(self._n_tile_rows): + for tc in range(self._n_tile_cols): + cx, cy = self._tile_centers[(tr, tc)] + dist = np.sqrt((cam_x - cx) ** 2 + (cam_y - cy) ** 2) + + lod = compute_lod_level(dist, self._lod_distances) + lod = min(lod, self._max_lod) + + tile_id = _tile_gid(tr, tc) + new_tile_ids.add(tile_id) + + prev_lod = self._tile_lods.get((tr, tc), -1) + if lod != prev_lod or force: + verts, indices = self._get_tile_mesh(tr, tc, lod) + if verts is not None: + # Apply VE + if ve != 1.0: + verts = verts.copy() + verts[2::3] *= ve + rtx.add_geometry(tile_id, verts, indices) + self._tile_lods[(tr, tc)] = lod + changed = True + + # Remove stale tiles (shouldn't happen, but be safe) + for old_id in self._active_tiles - new_tile_ids: + rtx.remove_geometry(old_id) + changed = True + + self._active_tiles = new_tile_ids + return changed + + def remove_all(self, rtx): + """Remove all LOD tile geometries from the scene.""" + for tile_id in list(self._active_tiles): + rtx.remove_geometry(tile_id) + self._active_tiles.clear() + self._tile_lods.clear() + self._last_update_pos = None + + def get_stats(self): + """Return a summary string of LOD state.""" + if not self._tile_lods: + return "LOD: no tiles" + from collections import Counter + counts = Counter(self._tile_lods.values()) + parts = [f"L{lvl}:{cnt}" for lvl, cnt in sorted(counts.items())] + total = self._n_tile_rows * self._n_tile_cols + return f"LOD tiles: {total} ({', '.join(parts)})" + + # ------------------------------------------------------------------ + # Internal + # ------------------------------------------------------------------ + + def _get_tile_mesh(self, tr, tc, lod): + """Build or retrieve cached tile mesh.""" + cache_key = (tr, tc, lod, self._base_subsample) + cached = self._tile_cache.get(cache_key) + if cached is not None: + return cached[0].copy(), cached[1].copy() + + verts, indices = self._build_tile_mesh(tr, tc, lod) + if verts is not None: + self._tile_cache[cache_key] = (verts.copy(), indices.copy()) + return verts, indices + + def _build_tile_mesh(self, tr, tc, lod): + """Triangulate a single tile at the given LOD level.""" + from .. import mesh as mesh_mod + + subsample = self._base_subsample * (2 ** lod) + r0 = tr * self._tile_size + c0 = tc * self._tile_size + r1 = min(r0 + self._tile_size, self._H) + c1 = min(c0 + self._tile_size, self._W) + + # Extract tile data with subsampling + tile = self._terrain_np[r0:r1:subsample, c0:c1:subsample] + th, tw = tile.shape + if th < 2 or tw < 2: + return None, None + + # Triangulate using the fast numba/CUDA path + num_verts = th * tw + num_tris = (th - 1) * (tw - 1) * 2 + verts = np.zeros(num_verts * 3, dtype=np.float32) + indices = np.zeros(num_tris * 3, dtype=np.int32) + mesh_mod.triangulate_terrain(verts, indices, tile, scale=1.0) + + # Transform from local grid coords to world coords. + # triangulate_terrain writes x=w, y=h in grid-local pixel indices. + # We need: x = (c0 + w*subsample) * psx + # y = (r0 + h*subsample) * psy + verts[0::3] = verts[0::3] * subsample * self._psx + c0 * self._psx + verts[1::3] = verts[1::3] * subsample * self._psy + r0 * self._psy + + # Add edge skirt to hide T-junction cracks between LOD levels + verts, indices = _add_tile_skirt(verts, indices, th, tw) + + return verts, indices + + +# ------------------------------------------------------------------ +# Module-level helpers +# ------------------------------------------------------------------ + +def _tile_gid(tr, tc): + """Geometry ID for a terrain LOD tile.""" + return f'terrain_lod_r{tr}_c{tc}' + + +def is_terrain_lod_gid(gid): + """Return True if *gid* belongs to a terrain LOD tile.""" + return gid.startswith('terrain_lod_r') + + +def _add_tile_skirt(vertices, indices, H, W, skirt_depth=None): + """Add a thin skirt around tile edges. + + The skirt is deliberately small — just enough to cover gaps at + LOD boundaries. It uses the same algorithm as + ``mesh.add_terrain_skirt`` but with a shallower default depth. + """ + z_vals = vertices[2::3] + z_min = float(np.nanmin(z_vals)) + z_max = float(np.nanmax(z_vals)) + + if skirt_depth is None: + z_range = z_max - z_min + skirt_depth = max(0.5, z_range * 0.02) + + skirt_z = z_min - skirt_depth + + # Build clockwise perimeter (same order as mesh.add_terrain_skirt) + top = np.arange(W, dtype=np.int32) + right = (np.arange(1, H, dtype=np.int32)) * W + (W - 1) + bottom = (H - 1) * W + np.arange(W - 2, -1, -1, dtype=np.int32) + left = np.arange(H - 2, 0, -1, dtype=np.int32) * W + perim = np.concatenate([top, right, bottom, left]) + n_perim = len(perim) + n_orig = len(vertices) // 3 + + skirt_verts = np.empty(n_perim * 3, dtype=np.float32) + skirt_verts[0::3] = vertices[perim * 3] + skirt_verts[1::3] = vertices[perim * 3 + 1] + skirt_verts[2::3] = skirt_z + + idx = np.arange(n_perim, dtype=np.int32) + idx_next = np.roll(idx, -1) + top_a = perim + top_b = perim[idx_next] + bot_a = (n_orig + idx).astype(np.int32) + bot_b = (n_orig + idx_next).astype(np.int32) + + wall_tris = np.empty(n_perim * 6, dtype=np.int32) + wall_tris[0::6] = top_a + wall_tris[1::6] = bot_b + wall_tris[2::6] = top_b + wall_tris[3::6] = top_a + wall_tris[4::6] = bot_a + wall_tris[5::6] = bot_b + + new_verts = np.concatenate([vertices, skirt_verts]) + new_indices = np.concatenate([indices, wall_tris]) + return new_verts, new_indices