fix(simulation): Layer-1 self-describing grid dims + Moore basin trace (#960)

Two Layer-1 generation fixes surfaced by the per-layer atlas viewer:

- Layer1Output now carries grid_w/grid_h (the downsampled working-grid the
  positions live in). The client maps overlays from these, so the scale is
  correct for any source heightmap resolution rather than assuming the texture
  size — fixes overlays projecting at half scale into a corner.
- Drainage basin boundaries are traced as ordered, non-self-crossing contours
  via Moore-neighbour tracing instead of an angle-from-centroid sort. The sort
  produced star-shaped, self-crossing polygons for concave basins that rendered
  as straight chords across the map.

Golden (cascade_layer1.json) and the cross-language atlas_response_ready
fixture regenerated. 100 atlas lib tests + the new tracer test pass.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

server/src/atlas/drainage.rs

@@ -543,43 +543,16 @@ fn build_basins(labels: &[i32], w: usize, h: usize) -> Vec<DrainageBasin> {
         let cells = &basin_map[&basin_id];
         let area_pct = cells.len() as f32 / n as f32;
 
-        // Boundary cells: in this basin, adjacent to a different basin or edge.
-        let mut boundary: Vec<(u16, u16)> = Vec::new();
-        for &idx in cells {
-            let r = idx / w;
-            let c = idx % w;
-            let mut on_boundary = false;
-            for &(dr, dc) in &D8 {
-                let nr = r as i32 + dr;
-                let nc = (c as i32 + dc).rem_euclid(w as i32) as usize;
-                if nr < 0 || nr >= h as i32 {
-                    on_boundary = true;
-                    break;
-                }
-                if labels[nr as usize * w + nc] != basin_id {
-                    on_boundary = true;
-                    break;
-                }
-            }
-            if on_boundary {
-                boundary.push((r as u16, c as u16));
-            }
-        }
-
-        // Sort boundary by angle from centroid for a coherent polygon.
-        if !boundary.is_empty() {
-            let cr = boundary.iter().map(|&(r, _)| r as f32).sum::<f32>() / boundary.len() as f32;
-            let cc = boundary.iter().map(|&(_, c)| c as f32).sum::<f32>() / boundary.len() as f32;
-            boundary.sort_by(|&(r1, c1), &(r2, c2)| {
-                let a1 = (r1 as f32 - cr).atan2(c1 as f32 - cc);
-                let a2 = (r2 as f32 - cr).atan2(c2 as f32 - cc);
-                a1.partial_cmp(&a2).unwrap_or(std::cmp::Ordering::Equal)
-            });
-            // Subsample to ≤500 points.
-            if boundary.len() > 500 {
-                let step = boundary.len() / 500;
-                boundary = boundary.into_iter().step_by(step).collect();
-            }
+        // Outer boundary as an ordered, non-self-crossing contour via Moore-
+        // neighbour tracing from the basin's first (row-major) cell. The previous
+        // angle-from-centroid sort produced star-shaped, self-crossing polygons for
+        // concave basins, which rendered as straight chords across the map (#960).
+        let start = cells.iter().copied().min().unwrap_or(0);
+        let mut boundary = trace_outer_boundary(labels, w, h, basin_id, start);
+        // Decimate to ≤500 points, preserving traversal order (and thus shape).
+        if boundary.len() > 500 {
+            let step = boundary.len() / 500;
+            boundary = boundary.into_iter().step_by(step).collect();
         }
 
         basins.push(DrainageBasin {
@@ -592,6 +565,67 @@ fn build_basins(labels: &[i32], w: usize, h: usize) -> Vec<DrainageBasin> {
     basins
 }
 
+/// Trace the outer boundary of the connected component of `basin_id` containing
+/// `start` (a row-major cell index), clockwise, via Moore-neighbour tracing.
+/// Produces an ordered, 8-connected, non-self-crossing perimeter. The grid edge is
+/// treated as background (no x-wrap) — this is for atlas visualization, not flow.
+fn trace_outer_boundary(labels: &[i32], w: usize, h: usize, basin_id: i32, start: usize) -> Vec<(u16, u16)> {
+    // Moore-neighbourhood offsets in clockwise order: N, NE, E, SE, S, SW, W, NW.
+    const DIRS: [(i32, i32); 8] = [
+        (-1, 0),
+        (-1, 1),
+        (0, 1),
+        (1, 1),
+        (1, 0),
+        (1, -1),
+        (0, -1),
+        (-1, -1),
+    ];
+    let is_fg = |r: i32, c: i32| -> bool {
+        r >= 0
+            && r < h as i32
+            && c >= 0
+            && c < w as i32
+            && labels[r as usize * w + c as usize] == basin_id
+    };
+    let dir_index = |dr: i32, dc: i32| -> usize { DIRS.iter().position(|&o| o == (dr, dc)).unwrap_or(0) };
+
+    let sr = (start / w) as i32;
+    let sc = (start % w) as i32;
+    let s = (sr, sc);
+    let mut boundary: Vec<(u16, u16)> = vec![(sr as u16, sc as u16)];
+    let mut p = s;
+    // Backtrack starts west of `start`: it is the first cell in scan order, so its
+    // western neighbour is background. Consecutive Moore neighbours are 8-adjacent,
+    // so the new backtrack stays adjacent to the new boundary cell each step.
+    let mut b = (sr, sc - 1);
+    let max_steps = w * h * 8 + 16;
+    for _ in 0..max_steps {
+        let b_idx = dir_index(b.0 - p.0, b.1 - p.1);
+        let mut prev = b;
+        let mut advanced = false;
+        for k in 1..=8 {
+            let d = (b_idx + k) % 8;
+            let cand = (p.0 + DIRS[d].0, p.1 + DIRS[d].1);
+            if is_fg(cand.0, cand.1) {
+                if cand == s {
+                    return boundary; // closed the loop (start already at index 0)
+                }
+                boundary.push((cand.0 as u16, cand.1 as u16));
+                b = prev;
+                p = cand;
+                advanced = true;
+                break;
+            }
+            prev = cand;
+        }
+        if !advanced {
+            break; // isolated cell — no foreground neighbour
+        }
+    }
+    boundary
+}
+
 // ---------------------------------------------------------------------------
 // Tests
 // ---------------------------------------------------------------------------
@@ -616,6 +650,33 @@ mod tests {
             .collect()
     }
 
+    #[test]
+    fn trace_outer_boundary_is_an_adjacent_contour() {
+        // Concave (L-shaped) basin (id=1) in a 4x4 grid; -1 is background.
+        let labels: Vec<i32> = vec![
+            1, 1, -1, -1, //
+            1, 1, -1, -1, //
+            1, 1, 1, 1, //
+            1, 1, 1, 1, //
+        ];
+        let boundary = trace_outer_boundary(&labels, 4, 4, 1, 0);
+        assert!(boundary.len() >= 8, "expected a real perimeter, got {boundary:?}");
+        // The defining property the angle-sort violated: consecutive boundary
+        // points are 8-adjacent (a genuine contour walk, not crossing chords).
+        for w in boundary.windows(2) {
+            let dr = (w[0].0 as i32 - w[1].0 as i32).abs();
+            let dc = (w[0].1 as i32 - w[1].1 as i32).abs();
+            assert!(
+                dr <= 1 && dc <= 1 && dr + dc > 0,
+                "non-adjacent step {:?} -> {:?}",
+                w[0],
+                w[1]
+            );
+        }
+        // Deterministic (D-010): same input → same trace.
+        assert_eq!(boundary, trace_outer_boundary(&labels, 4, 4, 1, 0));
+    }
+
     #[test]
     fn flat_grid_produces_single_basin() {
         let elev = flat_grid(16, 8, 0.5);

server/src/atlas/layer1.rs

@@ -32,6 +32,12 @@ pub struct Layer1Output {
     /// Geographic attractors (D-209) with sub-biome + cost (D-210), sorted by
     /// `(attractor_type, row, col)`.
     pub attractors: Vec<GeographicAttractor>,
+    /// Working-grid dimensions every position in this output (river cells, basin
+    /// boundaries, attractor positions) is expressed in — equals the downsampled
+    /// heightmap size. The client maps these onto the displayed heightmap (#960),
+    /// so the overlay scale stays correct for any source resolution (mod-safe).
+    pub grid_w: u32,
+    pub grid_h: u32,
 }
 
 /// Run the Layer-1 topography pipeline for a single body.
@@ -60,6 +66,8 @@ pub fn run_layer1(hm: &BodyHeightmap) -> Layer1Output {
         river_network: drainage.river_network,
         drainage_basins: drainage.drainage_basins,
         attractors,
+        grid_w: hm.width,
+        grid_h: hm.height,
     }
 }
 

server/src/atlas/layer_proxy.rs

@@ -70,6 +70,10 @@ pub fn handle_atlas_request(
             river_network: state.river_network.clone(),
             drainage_basins: state.drainage_basins.clone(),
             attractors: state.attractors.clone(),
+            // The cascade ran on the downsampled heightmap, so its dims are the
+            // working grid all Layer-1 positions are expressed in (#960).
+            grid_w: state.heightmap_width,
+            grid_h: state.heightmap_height,
         };
         return AtlasLayerResponse {
             body_id: req.body_id.clone(),

server/tests/gen_fixtures.rs

@@ -561,6 +561,8 @@ fn generate_atlas_layer_response_fixtures() {
             sub_biome: SubBiomeVariant::CoastalLowland,
             terrain_modification_cost: 1.7,
         }],
+        grid_w: 512,
+        grid_h: 256,
     };
     let ready = AtlasLayerResponse {
         body_id: "GJ1c".into(),