feat(world): render steep tiles on a smooth bicubic surface, sampled in parallel

Tiles holding a slope over 35 degrees are meshed at two vertices per cell
on a clamped Catmull-Rom surface that passes through every grid node, so
cliffs and ledges lose the facets of the heightfield grid while the
collision stays on the grid. The domain field exposes the smooth
elevation and a separable, per-phase-weighted block sampler with linear
edge extrapolation and a one-cell margin for continuous normals; the actor
samples tiles on the worker thread pool and builds meshes afterwards.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

docs/briefs/terrain-site-zone1-gameplay.md

@@ -245,7 +245,14 @@ de contacts existant les sépare du fond.
   à 45 m, indigo au pied du tombant).
 - Maillage découpé en **64 tuiles de 96 × 96 m** (64 × 64 cellules), avec
   `visibility_range` pour ne dessiner que les tuiles dans le brouillard
-  utile ; une seule `HeightMapShape3D` pour la collision.
+  utile ; une seule `HeightMapShape3D` pour la collision. Les tuiles qui
+  contiennent une pente de plus de 35° sont rendues à deux sommets par
+  cellule sur la surface bicubique (Catmull-Rom, qui passe par chaque nœud
+  du champ, bornée entre nœuds pour ne jamais dépasser une marche) : les
+  facettes de la grille disparaissent sur les parois, la collision reste sur
+  la grille et coïncide avec le rendu aux nœuds. L'échantillonnage des tuiles
+  est pur et s'exécute en parallèle (`WorkerThreadPool`), les maillages sont
+  assemblés ensuite sur le thread principal.
 - Props low-poly posés à la main, en nombre borné : blocs d'éboulis au pied
   du tombant, têtes de corail sur les patch reefs, rien d'animé.
 - Brouillard, lumière directionnelle et particules de la scène de test

game/domain/sites/terrain_height_field.gd

@@ -151,6 +151,193 @@ func elevation_at(x_m: float, z_m: float) -> float:
 	return lerpf(top, bottom, tz)
 
 
+## Smooth (Catmull-Rom bicubic) elevation: passes through every cell value, so it
+## agrees with elevation_at on the grid nodes and bends between them instead of
+## creasing. Reproduces planes exactly, including at the grid edges, where the field
+## is extrapolated linearly. Each pass is clamped between its two bracketing nodes,
+## so the surface never overshoots a step (no dip at the foot of a block, no bump
+## above a ledge). Used to render steep tiles at a finer step.
+func smooth_elevation_at(x_m: float, z_m: float) -> float:
+	var column_f := clampf((x_m - _data.origin_x_m) / _data.cell_size_m, 0.0, float(_data.columns - 1))
+	var row_f := clampf((z_m - _data.origin_z_m) / _data.cell_size_m, 0.0, float(_data.rows - 1))
+	var column := mini(int(floorf(column_f)), _data.columns - 2)
+	var row := mini(int(floorf(row_f)), _data.rows - 2)
+	var weights_x := _catmull_rom_weights(column_f - float(column))
+	var weights_z := _catmull_rom_weights(row_f - float(row))
+	var along := PackedFloat32Array()
+	along.resize(4)
+	for offset_z in 4:
+		var sample_row := row - 1 + offset_z
+		var p1 := _extrapolated_elevation(column, sample_row)
+		var p2 := _extrapolated_elevation(column + 1, sample_row)
+		var along_x := (
+			weights_x[0] * _extrapolated_elevation(column - 1, sample_row)
+			+ weights_x[1] * p1
+			+ weights_x[2] * p2
+			+ weights_x[3] * _extrapolated_elevation(column + 2, sample_row)
+		)
+		along[offset_z] = clampf(along_x, minf(p1, p2), maxf(p1, p2))
+	var result := weights_z[0] * along[0] + weights_z[1] * along[1] + weights_z[2] * along[2] + weights_z[3] * along[3]
+	return clampf(result, minf(along[1], along[2]), maxf(along[1], along[2]))
+
+
+## Positions and normals of one cell block on the smooth surface.
+class SurfaceBlock extends RefCounted:
+	var positions := PackedVector3Array()
+	var normals := PackedVector3Array()
+
+
+## The cell block [start_column..end_column] x [start_row..end_row] sampled on the smooth
+## surface every 1 / subdivision cell, row-major: (cells * subdivision + 1) points per
+## axis. The points on the grid nodes equal cell_world_position; the normals are
+## central differences on the smooth surface, continuous across tile borders thanks to
+## a one-cell margin. The smooth grid is built once for both arrays.
+func cell_block_subdivided(
+	start_column: int, end_column: int, start_row: int, end_row: int, subdivision: int
+) -> SurfaceBlock:
+	var block := SurfaceBlock.new()
+	if subdivision <= 1:
+		block.positions = cell_block_positions(start_column, end_column, start_row, end_row)
+		block.normals = cell_block_normals(start_column, end_column, start_row, end_row)
+		return block
+	var first_column := _clamp_column(start_column)
+	var last_column := _clamp_column(end_column)
+	var first_row := _clamp_row(start_row)
+	var last_row := _clamp_row(end_row)
+	var fine := _fine_block(first_column, last_column, first_row, last_row, subdivision)
+	var points_x := (last_column - first_column) * subdivision + 1
+	var points_z := (last_row - first_row) * subdivision + 1
+	var padded_x := points_x + 2 * subdivision
+	var step_m := _data.cell_size_m / float(subdivision)
+	var x0 := _data.origin_x_m + float(first_column) * _data.cell_size_m
+	var z0 := _data.origin_z_m + float(first_row) * _data.cell_size_m
+	var twice_step_m := 2.0 * step_m
+	var positions := PackedVector3Array()
+	var normals := PackedVector3Array()
+	positions.resize(points_x * points_z)
+	normals.resize(points_x * points_z)
+	var index := 0
+	for sub_row in points_z:
+		var fine_row_start := (sub_row + subdivision) * padded_x + subdivision
+		var z := z0 + float(sub_row) * step_m
+		for sub_column in points_x:
+			var at := fine_row_start + sub_column
+			positions[index] = Vector3(x0 + float(sub_column) * step_m, fine[at], z)
+			var slope_x := (fine[at + 1] - fine[at - 1]) / twice_step_m
+			var slope_z := (fine[at + padded_x] - fine[at - padded_x]) / twice_step_m
+			normals[index] = Vector3(-slope_x, 1.0, -slope_z).normalized()
+			index += 1
+	block.positions = positions
+	block.normals = normals
+	return block
+
+
+func cell_block_positions_subdivided(
+	start_column: int, end_column: int, start_row: int, end_row: int, subdivision: int
+) -> PackedVector3Array:
+	return cell_block_subdivided(start_column, end_column, start_row, end_row, subdivision).positions
+
+
+func cell_block_normals_subdivided(
+	start_column: int, end_column: int, start_row: int, end_row: int, subdivision: int
+) -> PackedVector3Array:
+	return cell_block_subdivided(start_column, end_column, start_row, end_row, subdivision).normals
+
+
+## Smooth heights of the block padded by one cell on every side, row-major with
+## (cells + 2) * subdivision + 1 points per axis. Separable Catmull-Rom with the
+## weights of each sub-cell phase computed once: no per-point function calls.
+func _fine_block(first_column: int, last_column: int, first_row: int, last_row: int, subdivision: int) -> PackedFloat32Array:
+	var cells_x := last_column - first_column + 2
+	var cells_z := last_row - first_row + 2
+	var source_x := cells_x + 4
+	var source_z := cells_z + 4
+	var source_column0 := first_column - 2
+	var source_row0 := first_row - 2
+	var fine_x := cells_x * subdivision + 1
+	var fine_z := cells_z * subdivision + 1
+	var phase_weights: Array[PackedFloat32Array] = []
+	for phase in subdivision:
+		phase_weights.append(_catmull_rom_weights(float(phase) / float(subdivision)))
+	# Source samples with linear extrapolation outside the grid.
+	var source := PackedFloat32Array()
+	source.resize(source_x * source_z)
+	for sz in source_z:
+		for sx in source_x:
+			source[sz * source_x + sx] = _extrapolated_elevation(source_column0 + sx, source_row0 + sz)
+	# Pass 1: along x, for every source row.
+	var along := PackedFloat32Array()
+	along.resize(fine_x * source_z)
+	for sz in source_z:
+		var row_start := sz * source_x
+		var out_start := sz * fine_x
+		for cell in cells_x:
+			var p0 := source[row_start + cell]
+			var p1 := source[row_start + cell + 1]
+			var p2 := source[row_start + cell + 2]
+			var p3 := source[row_start + cell + 3]
+			var low := minf(p1, p2)
+			var high := maxf(p1, p2)
+			for phase in subdivision:
+				var w := phase_weights[phase]
+				along[out_start + cell * subdivision + phase] = clampf(w[0] * p0 + w[1] * p1 + w[2] * p2 + w[3] * p3, low, high)
+		along[out_start + cells_x * subdivision] = source[row_start + cells_x + 1]
+	# Pass 2: along z, for every fine column.
+	var fine := PackedFloat32Array()
+	fine.resize(fine_x * fine_z)
+	for cell in cells_z:
+		for phase in subdivision:
+			var w := phase_weights[phase]
+			var out_row := (cell * subdivision + phase) * fine_x
+			var r0 := cell * fine_x
+			var r1 := r0 + fine_x
+			var r2 := r1 + fine_x
+			var r3 := r2 + fine_x
+			for fx in fine_x:
+				var a1 := along[r1 + fx]
+				var a2 := along[r2 + fx]
+				fine[out_row + fx] = clampf(
+					w[0] * along[r0 + fx] + w[1] * a1 + w[2] * a2 + w[3] * along[r3 + fx],
+					minf(a1, a2),
+					maxf(a1, a2)
+				)
+	var last_out := cells_z * subdivision * fine_x
+	var last_source := (cells_z + 1) * fine_x
+	for fx in fine_x:
+		fine[last_out + fx] = along[last_source + fx]
+	return fine
+
+
+## Cell elevation, extended past the grid by linear extrapolation so the smooth surface
+## keeps reproducing planes at the edges.
+func _extrapolated_elevation(column: int, row: int) -> float:
+	var columns := _data.columns
+	var rows := _data.rows
+	var c := clampi(column, 0, columns - 1)
+	var r := clampi(row, 0, rows - 1)
+	var value := _data.elevations_m[r * columns + c]
+	if column < 0:
+		value += float(column) * (_data.elevations_m[r * columns + 1] - _data.elevations_m[r * columns])
+	elif column > columns - 1:
+		value += float(column - columns + 1) * (_data.elevations_m[r * columns + columns - 1] - _data.elevations_m[r * columns + columns - 2])
+	if row < 0:
+		value += float(row) * (_data.elevations_m[columns + c] - _data.elevations_m[c])
+	elif row > rows - 1:
+		value += float(row - rows + 1) * (_data.elevations_m[(rows - 1) * columns + c] - _data.elevations_m[(rows - 2) * columns + c])
+	return value
+
+
+static func _catmull_rom_weights(t: float) -> PackedFloat32Array:
+	var t2 := t * t
+	var t3 := t2 * t
+	return PackedFloat32Array([
+		0.5 * (-t + 2.0 * t2 - t3),
+		0.5 * (2.0 - 5.0 * t2 + 3.0 * t3),
+		0.5 * (t + 4.0 * t2 - 3.0 * t3),
+		0.5 * (-t2 + t3),
+	])
+
+
 func material_at(x_m: float, z_m: float) -> int:
 	var column := int(roundf((x_m - _data.origin_x_m) / _data.cell_size_m))
 	var row := int(roundf((z_m - _data.origin_z_m) / _data.cell_size_m))

game/world/sites/bathymetric_terrain_actor.gd

@@ -8,6 +8,11 @@ const TERRAIN_SHADER = preload("res://game/world/sites/seabed_terrain.gdshader")
 const TILE_CELLS := 64
 const VISIBILITY_RANGE_END_M := 420.0
 const VISIBILITY_RANGE_MARGIN_M := 60.0
+## Tiles holding a cell steeper than this are rendered on the smooth surface at a finer
+## step, which removes the facets of the heightfield grid on cliffs and ledges. The
+## collision stays on the grid: the smooth surface passes through every grid node.
+const STEEP_TILE_SLOPE_DEGREES := 35.0
+const STEEP_TILE_SUBDIVISION := 2
 const DETAIL_TEXTURE_SEED := 4102
 
 var _field: TerrainHeightField
@@ -46,6 +51,14 @@ func tile_count() -> int:
 	return 0 if _tiles_root == null else _tiles_root.get_child_count()
 
 
+## Mesh points per cell edge of a tile: 1 for a flat tile, STEEP_TILE_SUBDIVISION for a
+## steep one.
+func tile_subdivision(tile_index: int) -> int:
+	if _tiles_root == null or tile_index < 0 or tile_index >= _tiles_root.get_child_count():
+		return 0
+	return int(_tiles_root.get_child(tile_index).get_meta(&"subdivision", 1))
+
+
 func sonar_material_at(world_position: Vector3) -> int:
 	return _field.material_at(world_position.x, world_position.z)
 
@@ -65,10 +78,20 @@ func _clear() -> void:
 		_collision = null
 
 
+## The geometry of one tile, computed off the main thread, then turned into a mesh.
+class TileGeometry extends RefCounted:
+	var cells := Rect2i()  # position = first column / row, end = last column / row
+	var subdivision := 1
+	var centre := Vector3.ZERO
+	var vertices := PackedVector3Array()
+	var normals := PackedVector3Array()
+
+
 func _build_tiles() -> void:
 	_tiles_root = Node3D.new()
 	_tiles_root.name = "Tiles"
 	add_child(_tiles_root)
+	var tiles: Array[TileGeometry] = []
 	var last_column := _field.columns() - 1
 	var last_row := _field.rows() - 1
 	var start_row := 0
@@ -77,27 +100,55 @@ func _build_tiles() -> void:
 		var start_column := 0
 		while start_column < last_column:
 			var end_column := mini(start_column + TILE_CELLS, last_column)
-			_tiles_root.add_child(_build_tile(start_column, end_column, start_row, end_row))
+			var tile := TileGeometry.new()
+			tile.cells = Rect2i(start_column, start_row, end_column - start_column, end_row - start_row)
+			tiles.append(tile)
 			start_column = end_column
 		start_row = end_row
+	# The per-tile sampling is pure (reads the field only), so the tiles are computed in
+	# parallel; meshes and nodes are then built on the main thread in tile order.
+	var group := WorkerThreadPool.add_group_task(
+		_compute_tile_geometry.bind(tiles), tiles.size(), -1, true, "seabed tiles"
+	)
+	WorkerThreadPool.wait_for_group_task_completion(group)
+	for tile in tiles:
+		_tiles_root.add_child(_build_tile(tile))
 
 
-func _build_tile(start_column: int, end_column: int, start_row: int, end_row: int) -> MeshInstance3D:
+func _compute_tile_geometry(tile_index: int, tiles: Array[TileGeometry]) -> void:
+	var tile := tiles[tile_index]
+	var start_column := tile.cells.position.x
+	var start_row := tile.cells.position.y
+	var end_column := tile.cells.end.x
+	var end_row := tile.cells.end.y
 	var corner_min := _field.cell_world_position(start_column, start_row)
 	var corner_max := _field.cell_world_position(end_column, end_row)
-	var centre := Vector3(
+	tile.centre = Vector3(
 		(corner_min.x + corner_max.x) * 0.5,
 		0.0,
 		(corner_min.z + corner_max.z) * 0.5
 	)
-	var vertices := _field.cell_block_positions(start_column, end_column, start_row, end_row)
-	for index in vertices.size():
-		vertices[index] -= centre
 	var normals := _field.cell_block_normals(start_column, end_column, start_row, end_row)
+	tile.subdivision = STEEP_TILE_SUBDIVISION if _is_steep(normals) else 1
+	var vertices: PackedVector3Array
+	if tile.subdivision > 1:
+		var block := _field.cell_block_subdivided(start_column, end_column, start_row, end_row, tile.subdivision)
+		vertices = block.positions
+		normals = block.normals
+	else:
+		vertices = _field.cell_block_positions(start_column, end_column, start_row, end_row)
+	for index in vertices.size():
+		vertices[index] -= tile.centre
+	tile.vertices = vertices
+	tile.normals = normals
+
+
+func _build_tile(tile: TileGeometry) -> MeshInstance3D:
+	var subdivision := tile.subdivision
 	var indices := PackedInt32Array()
-	var columns_in_tile := end_column - start_column + 1
-	for row in range(end_row - start_row):
-		for column in range(end_column - start_column):
+	var columns_in_tile := tile.cells.size.x * subdivision + 1
+	for row in range(tile.cells.size.y * subdivision):
+		for column in range(tile.cells.size.x * subdivision):
 			var top_left := row * columns_in_tile + column
 			var top_right := top_left + 1
 			var bottom_left := top_left + columns_in_tile
@@ -106,22 +157,31 @@ func _build_tile(start_column: int, end_column: int, start_row: int, end_row: in
 			indices.append_array([top_right, bottom_right, bottom_left])
 	var arrays: Array = []
 	arrays.resize(Mesh.ARRAY_MAX)
-	arrays[Mesh.ARRAY_VERTEX] = vertices
-	arrays[Mesh.ARRAY_NORMAL] = normals
+	arrays[Mesh.ARRAY_VERTEX] = tile.vertices
+	arrays[Mesh.ARRAY_NORMAL] = tile.normals
 	arrays[Mesh.ARRAY_INDEX] = indices
 	var mesh := ArrayMesh.new()
 	mesh.add_surface_from_arrays(Mesh.PRIMITIVE_TRIANGLES, arrays)
 	var instance := MeshInstance3D.new()
-	instance.name = "Tile_%d_%d" % [start_column / TILE_CELLS, start_row / TILE_CELLS]
+	instance.name = "Tile_%d_%d" % [tile.cells.position.x / TILE_CELLS, tile.cells.position.y / TILE_CELLS]
 	instance.mesh = mesh
+	instance.set_meta(&"subdivision", subdivision)
 	instance.material_override = _material
-	instance.position = centre
+	instance.position = tile.centre
 	instance.visibility_range_end = VISIBILITY_RANGE_END_M
 	instance.visibility_range_end_margin = VISIBILITY_RANGE_MARGIN_M
 	instance.visibility_range_fade_mode = GeometryInstance3D.VISIBILITY_RANGE_FADE_SELF
 	return instance
 
 
+static func _is_steep(normals: PackedVector3Array) -> bool:
+	var steep_cosine := cos(deg_to_rad(STEEP_TILE_SLOPE_DEGREES))
+	for normal in normals:
+		if normal.y < steep_cosine:
+			return true
+	return false
+
+
 func _build_collision() -> void:
 	var shape := HeightMapShape3D.new()
 	shape.map_width = _field.columns()

tests/integration/run_site_terrain_integration.gd

@@ -19,7 +19,7 @@ const DIVE_SITE_SCENE = preload("res://game/world/sites/dive_site.tscn")
 const ZONE1_SITE_PATH := "res://content/sites/zone1_vigilant.tres"
 const SHOWCASE_SITE_PATH := "res://content/sites/terrain_showcase.tres"
 const GROOVE_PROBE_OFFSET_M := 25.0
-const TEST_COUNT := 41
+const TEST_COUNT := 44
 
 var _failures: Array[String] = []
 
@@ -109,12 +109,21 @@ func _test_actor_matches_field() -> void:
 		tile_vertices[tile_indices[2]]
 	)
 	_expect(first_face.normal.y > 0.0, "tile triangles are wound clockwise from above (front face up)")
-	# The first tile spans cells 0..64 on both axes, so vertex index = row * 65 + column.
+	# The first tile holds the steep bump, so it is rendered subdivided; the second tile
+	# (cells 64..80) is flat and stays on the grid.
+	var subdivision := actor.tile_subdivision(0)
+	_expect(subdivision == BathymetricTerrainActor.STEEP_TILE_SUBDIVISION, "the steep tile is subdivided (%d)" % subdivision)
+	_expect(actor.tile_subdivision(1) == 1, "the flat tile keeps one vertex per cell")
+	# The first tile spans cells 0..64 on both axes: vertex index = (row * n) * (64 n + 1) + column * n.
+	var points_per_row := 64 * subdivision + 1
 	var worst_vertex_error := 0.0
 	for cell: Vector2i in [Vector2i(0, 0), Vector2i(10, 50), Vector2i(64, 64)]:
-		var vertex := first_tile.global_transform * tile_vertices[cell.y * 65 + cell.x]
+		var vertex := first_tile.global_transform * tile_vertices[(cell.y * subdivision) * points_per_row + cell.x * subdivision]
 		worst_vertex_error = maxf(worst_vertex_error, vertex.distance_to(field.cell_world_position(cell.x, cell.y)))
-	_expect(worst_vertex_error < 0.01, "tile vertices sit on their field cells (worst %.4f m)" % worst_vertex_error)
+	_expect(worst_vertex_error < 0.01, "tile vertices on grid nodes sit on their field cells (worst %.4f m)" % worst_vertex_error)
+	var mid_vertex := first_tile.global_transform * tile_vertices[(50 * subdivision + 1) * points_per_row + 10 * subdivision + 1]
+	var mid_expected := field.smooth_elevation_at(mid_vertex.x, mid_vertex.z)
+	_expect(absf(mid_vertex.y - mid_expected) < 0.01, "subdivided vertices sit on the smooth surface (%.3f vs %.3f)" % [mid_vertex.y, mid_expected])
 	var space := actor.get_world_3d().direct_space_state
 	var worst_error := 0.0
 	for cell: Vector2i in [

tests/unit/run_site_domain_tests.gd

@@ -7,7 +7,7 @@ const SitePlacement = preload("res://game/domain/sites/site_placement.gd")
 const SiteRoute = preload("res://game/domain/sites/site_route.gd")
 const DiveSiteDefinition = preload("res://game/domain/sites/dive_site_definition.gd")
 
-const TEST_COUNT := 41
+const TEST_COUNT := 46
 
 var _failures: Array[String] = []
 
@@ -97,6 +97,38 @@ func _test_height_field_sampling() -> void:
 	_expect(field.material_at(-10.4, 100.0) == SonarMaterial.Value.ROCK, "material is the nearest cell")
 	_expect(field.material_at(-6.0, 103.0) == SonarMaterial.Value.SAND, "material of an interior cell")
 	_expect(field.normal_at_cell(1, 1).is_equal_approx(Vector3(0.0, 1.0, 1.0).normalized()), "normal leans against the downhill direction")
+	_expect(is_equal_approx(field.smooth_elevation_at(-10.0, 100.0), field.elevation_at(-10.0, 100.0)), "smooth elevation passes through the grid nodes")
+	var plane_data := SiteTerrainData.new()
+	plane_data.columns = 7
+	plane_data.rows = 6
+	plane_data.cell_size_m = 2.0
+	plane_data.origin_x_m = 0.0
+	plane_data.origin_z_m = 0.0
+	var plane_elevations := PackedFloat32Array()
+	var plane_materials := PackedByteArray()
+	for row in 6:
+		for column in 7:
+			plane_elevations.append(-20.0 + 0.5 * float(column) * 2.0 - 0.25 * float(row) * 2.0)
+			plane_materials.append(SonarMaterial.Value.SAND)
+	plane_data.elevations_m = plane_elevations
+	plane_data.materials = plane_materials
+	var plane := TerrainHeightField.from_data(plane_data)
+	_expect(
+		is_equal_approx(plane.smooth_elevation_at(5.0, 3.0), -20.0 + 2.5 - 0.75)
+		and is_equal_approx(plane.smooth_elevation_at(0.5, 0.5), -20.0 + 0.25 - 0.125),
+		"smooth elevation reproduces a plane between nodes, edges included"
+	)
+	var fine := field.cell_block_positions_subdivided(1, 2, 0, 3, 2)
+	var fine_normals := field.cell_block_normals_subdivided(1, 2, 0, 3, 2)
+	_expect(fine.size() == 3 * 7 and fine_normals.size() == 3 * 7, "a subdivided block has cells * 2 + 1 points per axis")
+	var nodes_match := fine.size() == 21
+	for row in 4:
+		for column in range(1, 3):
+			var index := (row * 2) * 3 + (column - 1) * 2
+			if nodes_match and not fine[index].is_equal_approx(field.cell_world_position(column, row)):
+				nodes_match = false
+	_expect(nodes_match, "subdivided block points on the grid nodes equal the cell positions")
+	_expect(fine_normals[0].is_normalized() and fine_normals[10].y > 0.0, "subdivided normals are unit length and point up")
 	var block_positions := field.cell_block_positions(1, 2, 0, 3)
 	var block_normals := field.cell_block_normals(1, 2, 0, 3)
 	var block_matches := block_positions.size() == 8 and block_normals.size() == 8