Terrain LOD refactor: edge stitching, hybrid mesh LOD, performance fixes (#87)

* Terrain LOD refactor: edge stitching, hybrid mesh LOD, performance fixes

Major enhancements to the terrain LOD system:

- Replace tile skirts with boundary vertex stitching (_stitch_tile_boundary)
  that interpolates Z from coarser neighbor's pyramid data, eliminating
  T-junction cracks without extra geometry
- Hybrid LOD for placed geometry: chunk manager applies simplify_mesh()
  during merge based on chunk LOD level (LOD 0=full, LOD 3=10% triangles)
- Extract update() into focused helper phases (A-E) reducing main method
  from 324 to ~194 lines
- Fix critical NameError: add missing cache_key in _rebuild_vertical_exaggeration
- Add build retry budget (max 3 attempts) with logging for failed tile builds
- Rate-limit I/O prefetches (max 8/tick) to prevent thread pool flooding
- Early frustum culling before LOD/roughness computation (~40-60% fewer
  evaluations for tiles behind camera)
- Single-pass future harvesting in _collect_completed_builds
- Pyramid cache eviction for levels above max_lod
- In-place simplify cache eviction (del loop vs dict rebuild)
- O(1) stale tile eviction via set lookup
- Simplified Z re-snapping: final_z = stored_z * ve
- GPU re-snap threshold bumped from 1K to 10K verts
- Float32 precision for compute_tile_roughness bilinear interpolation
- Remove dead _add_tile_skirt function and references
- Add get_metrics() structured API for programmatic LOD monitoring
- Rename _has_pending to _has_in_flight_work for clarity
- 114 tests (up from ~30), including edge stitching Z-value verification,
  hysteresis dead zone coverage, retry budget, and mesh simplification

* general sync up

* Fix ocean-fill NaN replacement crashing on dask-backed rasters

Use the already-computed numpy array (base_np) instead of copying
the raw raster data, which may be a dask array that doesn't support
boolean mask assignment.

* quick fix for scenes

Author: brendancol

.gitignore

@@ -45,3 +45,4 @@ examples/.ipynb_checkpoints/
 conda-output/
 examples/cache/
 *gtfs.json
+.claude/worktrees/

cuda/common.h

@@ -68,4 +68,6 @@ struct Params
     int                    _pad0;              // padding for pointer alignment
     // --- point cloud fields (offset 88) ---
     float*                 point_colors;       // per-point RGBA (4 floats per point, indexed by primitive_id)
+    // --- smooth normal fields (offset 96) ---
+    unsigned long long*    smooth_normal_table; // [2*instanceId]=normals_ptr, [2*instanceId+1]=indices_ptr
 };

cuda/kernel.cu

@@ -99,17 +99,46 @@ extern "C" __global__ void __closesthit__chit()
 
     float3 n;
     if (optixIsTriangleHit()) {
-        float3 data[3];
-        // Always use the 4-parameter overload for backward compatibility.
-        // The parameterless overload (OptiX 9.1+) requires ABI version 99,
-        // which needs driver 570+.  The 4-param form works on all versions.
-        OptixTraversableHandle gas = optixGetGASTraversableHandle();
-        unsigned int sbtIdx = optixGetSbtGASIndex();
-        float time = optixGetRayTime();
-        optixGetTriangleVertexData(gas, primIdx, sbtIdx, time, data);
-        float3 AB = data[1] - data[0];
-        float3 AC = data[2] - data[0];
-        n = normalize(cross(AB, AC));
+        // Check for per-vertex smooth normals for this instance
+        bool has_smooth = false;
+        if (params.smooth_normal_table != 0) {
+            unsigned int instId = optixGetInstanceId();
+            unsigned long long normals_ptr = params.smooth_normal_table[2 * instId];
+            if (normals_ptr != 0) {
+                has_smooth = true;
+                unsigned long long indices_ptr = params.smooth_normal_table[2 * instId + 1];
+                const float* norms = reinterpret_cast<const float*>(normals_ptr);
+                const int*   idx   = reinterpret_cast<const int*>(indices_ptr);
+
+                // Barycentric interpolation of vertex normals
+                const float2 bary = optixGetTriangleBarycentrics();
+                const float w0 = 1.0f - bary.x - bary.y;
+                const float w1 = bary.x;
+                const float w2 = bary.y;
+
+                const int i0 = idx[primIdx * 3];
+                const int i1 = idx[primIdx * 3 + 1];
+                const int i2 = idx[primIdx * 3 + 2];
+
+                n = normalize(make_float3(
+                    w0 * norms[i0 * 3]     + w1 * norms[i1 * 3]     + w2 * norms[i2 * 3],
+                    w0 * norms[i0 * 3 + 1] + w1 * norms[i1 * 3 + 1] + w2 * norms[i2 * 3 + 1],
+                    w0 * norms[i0 * 3 + 2] + w1 * norms[i1 * 3 + 2] + w2 * norms[i2 * 3 + 2]
+                ));
+            }
+        }
+
+        if (!has_smooth) {
+            // Flat shading fallback: face normal from triangle vertices
+            float3 data[3];
+            OptixTraversableHandle gas = optixGetGASTraversableHandle();
+            unsigned int sbtIdx = optixGetSbtGASIndex();
+            float time = optixGetRayTime();
+            optixGetTriangleVertexData(gas, primIdx, sbtIdx, time, data);
+            float3 AB = data[1] - data[0];
+            float3 AC = data[2] - data[0];
+            n = normalize(cross(AB, AC));
+        }
     } else {
         // Round curve tube: use face-up normal for terrain roads/rivers
         n = make_float3(0.0f, 0.0f, 1.0f);