feat(tools): irregular strata and angular boulder blocks

Strata gain a dip, a slow seeded warp that bends and notches the beds, and
a per-bed thickness and hardness jitter with ordered boundaries, so walls
never invert. BoulderField gains an angularity parameter: elongated
super-ellipse slabs with flat tops, sharp shoulders, random heading and
tilt, footprint area preserved. Both sites are re-baked with the new
settings; the coral heads stay rounded.

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

tools/site_terrain_baker.gd

@@ -309,14 +309,31 @@ class Outcrop extends Primitive:
 		return NO_MATERIAL
 
 
-## count half-ellipsoid boulders scattered uniformly in a disc by a seeded generator.
-## Boulders combine by maximum, so overlapping blocks never pile into one mound, and every
-## cell they cover is rock.
+## count boulders scattered uniformly in a disc by a seeded generator. Boulders combine
+## by maximum, so overlapping blocks never pile into one mound, and every cell they cover
+## is rock. angularity 0 gives half-ellipsoids; towards 1 the blocks become elongated
+## super-ellipses with flat tops, sharp shoulders, a random heading and a tilt, like
+## fallen slabs. The footprint area is preserved whatever the elongation.
 class BoulderField extends Primitive:
+	const MAXIMUM_ASPECT := 2.2
+	const MAXIMUM_TILT := 0.35
+
+	## One block: centre, radius, heading, elongation, outline exponent, shoulder
+	## exponent and tilt across its long axis.
+	class Block extends RefCounted:
+		var centre: Vector2
+		var radius_m: float
+		var yaw: float
+		var aspect: float
+		var outline_power: float
+		var shoulder_power: float
+		var tilt: float
+
 	var centre: Vector2
 	var radius_m: float
 	var height_ratio: float
-	var _boulders := PackedVector3Array()  # x, z, radius
+	var angularity: float
+	var _blocks: Array[Block] = []
 	var _reach_m: float
 
 	func _init(
@@ -326,40 +343,69 @@ class BoulderField extends Primitive:
 		p_count: int,
 		p_boulder_radius_min_m: float,
 		p_boulder_radius_max_m: float,
-		p_height_ratio: float
+		p_height_ratio: float,
+		p_angularity: float = 0.0
 	) -> void:
 		centre = p_centre
 		radius_m = p_radius_m
 		height_ratio = p_height_ratio
-		_reach_m = p_radius_m + p_boulder_radius_max_m
+		angularity = clampf(p_angularity, 0.0, 1.0)
+		_reach_m = p_radius_m + p_boulder_radius_max_m * sqrt(lerpf(1.0, MAXIMUM_ASPECT, angularity))
 		var generator := RandomNumberGenerator.new()
 		generator.seed = p_seed
 		for _index in p_count:
 			var distance := p_radius_m * sqrt(generator.randf())
 			var angle := TAU * generator.randf()
-			var boulder_radius := generator.randf_range(p_boulder_radius_min_m, p_boulder_radius_max_m)
-			_boulders.append(Vector3(
-				p_centre.x + distance * cos(angle),
-				p_centre.y + distance * sin(angle),
-				boulder_radius
-			))
+			var block := Block.new()
+			block.centre = p_centre + Vector2(distance * cos(angle), distance * sin(angle))
+			block.radius_m = generator.randf_range(p_boulder_radius_min_m, p_boulder_radius_max_m)
+			# Draw every shape number even at angularity 0 so a field keeps the same
+			# positions and radii when its angularity is changed.
+			var yaw := TAU * generator.randf()
+			var aspect := generator.randf_range(1.0, MAXIMUM_ASPECT)
+			var outline := generator.randf_range(3.0, 6.0)
+			var shoulder := generator.randf_range(3.0, 6.0)
+			var tilt := generator.randf_range(-MAXIMUM_TILT, MAXIMUM_TILT)
+			block.yaw = yaw
+			block.aspect = lerpf(1.0, aspect, angularity)
+			block.outline_power = lerpf(2.0, outline, angularity)
+			block.shoulder_power = lerpf(2.0, shoulder, angularity)
+			block.tilt = tilt * angularity
+			_blocks.append(block)
 
+	## Centre x, centre z and radius of every block, in generation order.
 	func boulders() -> PackedVector3Array:
-		return _boulders
+		var result := PackedVector3Array()
+		for block in _blocks:
+			result.append(Vector3(block.centre.x, block.centre.y, block.radius_m))
+		return result
 
 	func elevation_delta(x_m: float, z_m: float) -> float:
 		if Vector2(x_m, z_m).distance_to(centre) > _reach_m:
 			return 0.0
 		var highest := 0.0
-		for boulder in _boulders:
-			var dx := x_m - boulder.x
-			var dz := z_m - boulder.y
-			var squared := dx * dx + dz * dz
-			var radius_squared := boulder.z * boulder.z
-			if squared < radius_squared:
-				highest = maxf(highest, boulder.z * height_ratio * sqrt(1.0 - squared / radius_squared))
+		for block in _blocks:
+			highest = maxf(highest, _block_height(block, x_m, z_m))
 		return highest
 
+	func _block_height(block: Block, x_m: float, z_m: float) -> float:
+		var dx := x_m - block.centre.x
+		var dz := z_m - block.centre.y
+		var long_axis_m := block.radius_m * sqrt(block.aspect)
+		if dx * dx + dz * dz >= long_axis_m * long_axis_m:
+			return 0.0
+		var short_axis_m := block.radius_m / sqrt(block.aspect)
+		var along := (cos(block.yaw) * dx + sin(block.yaw) * dz) / long_axis_m
+		var across := (-sin(block.yaw) * dx + cos(block.yaw) * dz) / short_axis_m
+		var outline := pow(
+			pow(absf(along), block.outline_power) + pow(absf(across), block.outline_power),
+			1.0 / block.outline_power
+		)
+		if outline >= 1.0:
+			return 0.0
+		var profile := pow(1.0 - pow(outline, block.shoulder_power), 1.0 / block.shoulder_power)
+		return maxf(block.radius_m * height_ratio * profile * (1.0 + block.tilt * along), 0.0)
+
 	func material_at(x_m: float, z_m: float) -> int:
 		if elevation_delta(x_m, z_m) > 0.0:
 			return SonarMaterial.Value.ROCK
@@ -444,32 +490,67 @@ class Modifier extends RefCounted:
 		return elevation_m
 
 
-## Quantises steep ground into horizontal ledges layer_height_m apart: the stratified
-## walls of a canyon or cliff. Flat floors and summits are below minimum_slope_degrees and
-## stay untouched; the effect is full from full_slope_degrees. sharpness 0 leaves the
-## field unchanged, 1 cuts hard steps. region_m is a world-space (x, z) rectangle with a
-## hard edge; an empty rectangle applies everywhere. Ledges need two to three cells of
-## wall per layer to be resolved: at 1.5 m per cell, 4 m layers on a 50 degree wall.
+## Quantises steep ground into ledges about layer_height_m apart: the stratified walls
+## of a canyon or cliff. Flat floors and summits are below minimum_slope_degrees and stay
+## untouched; the effect is full from full_slope_degrees. sharpness 0 leaves the field
+## unchanged, 1 cuts hard steps. region_m is a world-space (x, z) rectangle with a hard
+## edge; an empty rectangle applies everywhere.
+##
+## Real beds are neither horizontal nor equal: dip tilts the whole stack (metres of
+## elevation per metre along x and z), warp_amplitude_m bends the beds with slow seeded
+## noise so ledges pinch, swell and notch, and thickness_variation (0 to under 1) gives
+## every bed its own thickness and hardness, constant along the bed. Bed boundaries stay
+## ordered, so the quantised wall never inverts. All three default to the regular case.
+## Ledges need two to three cells of wall per layer to be resolved: at 1.5 m per cell,
+## 4 m layers on a 50 degree wall.
 class Strata extends Modifier:
+	const BED_SHARPNESS_SPREAD := 0.3
+
 	var layer_height_m: float
 	var sharpness: float
 	var region_m: Rect2
+	var dip: Vector2
+	var warp_amplitude_m: float
+	var thickness_variation: float
 	var _minimum_tangent: float
 	var _full_tangent: float
+	var _slow_warp: FastNoiseLite
+	var _notch_warp: FastNoiseLite
+	var _bed_noise: FastNoiseLite
 
 	func _init(
 		p_layer_height_m: float,
 		p_sharpness: float,
 		p_minimum_slope_degrees: float,
 		p_full_slope_degrees: float,
-		p_region_m: Rect2
+		p_region_m: Rect2,
+		p_dip: Vector2 = Vector2.ZERO,
+		p_warp_seed: int = 0,
+		p_warp_amplitude_m: float = 0.0,
+		p_thickness_variation: float = 0.0
 	) -> void:
 		assert(p_layer_height_m > 0.0, "strata need a positive layer height")
 		layer_height_m = p_layer_height_m
 		sharpness = clampf(p_sharpness, 0.0, 0.999)
 		region_m = p_region_m
+		dip = p_dip
+		warp_amplitude_m = p_warp_amplitude_m
+		thickness_variation = clampf(p_thickness_variation, 0.0, 0.9)
 		_minimum_tangent = tan(deg_to_rad(p_minimum_slope_degrees))
 		_full_tangent = tan(deg_to_rad(p_full_slope_degrees))
+		if p_warp_amplitude_m > 0.0:
+			_slow_warp = _make_noise(p_warp_seed, 0.016, 2)
+			_notch_warp = _make_noise(p_warp_seed + 1, 0.09, 2)
+		if thickness_variation > 0.0:
+			_bed_noise = _make_noise(p_warp_seed + 2, 1.0, 1)
+
+	static func _make_noise(seed: int, frequency: float, octaves: int) -> FastNoiseLite:
+		var noise := FastNoiseLite.new()
+		noise.noise_type = FastNoiseLite.TYPE_SIMPLEX_SMOOTH
+		noise.seed = seed
+		noise.frequency = frequency
+		noise.fractal_octaves = octaves
+		return noise
 
 	func modify(x_m: float, z_m: float, elevation_m: float, gradient: Vector2) -> float:
 		if region_m.has_area() and not region_m.has_point(Vector2(x_m, z_m)):
@@ -477,10 +558,43 @@ class Strata extends Modifier:
 		var weight := smoothstep(_minimum_tangent, _full_tangent, gradient.length())
 		if weight <= 0.0:
 			return elevation_m
-		var layer := floorf(elevation_m / layer_height_m)
-		var fraction := elevation_m / layer_height_m - layer
-		var sharpened := clampf((fraction - 0.5) / (1.0 - sharpness) + 0.5, 0.0, 1.0)
-		return lerpf(elevation_m, (layer + sharpened) * layer_height_m, weight)
+		# Work in "bed space": the stack tilted by the dip and bent by the warp.
+		var bed_elevation := elevation_m - dip.x * x_m - dip.y * z_m
+		if _slow_warp != null:
+			bed_elevation += _slow_warp.get_noise_2d(x_m, z_m) * warp_amplitude_m
+			bed_elevation += _notch_warp.get_noise_2d(x_m, z_m) * warp_amplitude_m * 0.25
+		var bed := floorf(bed_elevation / layer_height_m)
+		var floor_m := _bed_boundary(bed)
+		var ceiling_m := _bed_boundary(bed + 1.0)
+		# Jittered boundaries can leave the point just outside its nominal bed.
+		if bed_elevation < floor_m:
+			bed -= 1.0
+			ceiling_m = floor_m
+			floor_m = _bed_boundary(bed)
+		elif bed_elevation >= ceiling_m:
+			bed += 1.0
+			floor_m = ceiling_m
+			ceiling_m = _bed_boundary(bed + 1.0)
+		var fraction := (bed_elevation - floor_m) / (ceiling_m - floor_m)
+		var bed_sharpness := sharpness
+		if _bed_noise != null:
+			bed_sharpness = clampf(
+				sharpness + _bed_noise.get_noise_1d(bed * 17.0 + 5.0) * BED_SHARPNESS_SPREAD,
+				0.0,
+				0.999
+			)
+		var sharpened := clampf((fraction - 0.5) / (1.0 - bed_sharpness) + 0.5, 0.0, 1.0)
+		var quantised := floor_m + sharpened * (ceiling_m - floor_m)
+		# Apply the change in bed space back to the real elevation.
+		return elevation_m + (quantised - bed_elevation) * weight
+
+	## Lower boundary of a bed in bed space. The jitter is under half a layer, so the
+	## boundaries keep their order and every bed keeps a positive thickness.
+	func _bed_boundary(bed: float) -> float:
+		var boundary := bed * layer_height_m
+		if _bed_noise != null:
+			boundary += _bed_noise.get_noise_1d(bed * 17.0) * 0.5 * thickness_variation * layer_height_m
+		return boundary
 
 
 ## Seeded noise weighted by the local slope: rock faces and ridge flanks get crumbly

tools/sites/terrain_showcase_recipe.gd

@@ -68,16 +68,17 @@ static func recipe() -> SiteTerrainBaker.Recipe:
 		SiteTerrainBaker.Outcrop.new(EAST_BUTTE_CENTRE, 38.0, 75.0, 9.0, 8.0, 202, 5.0, true),
 		SiteTerrainBaker.Canyon.new(PackedVector2Array(CREVASSE_PATH), CREVASSE_HALF_WIDTH_M, CREVASSE_DEPTH_M, 1.0, 203, 1.0, SonarMaterial.Value.SAND),
 
-		# SW — le Chaos: ninety boulders of 2-6 m radius plus six coral heads, on rippled sand.
-		SiteTerrainBaker.BoulderField.new(BOULDER_FIELD_CENTRE, BOULDER_FIELD_RADIUS_M, 301, 90, 2.0, 6.0, 0.7),
-		SiteTerrainBaker.BoulderField.new(Vector2(-100.0, 130.0), 40.0, 302, 6, 6.0, 9.0, 0.6),
+		# SW — le Chaos: ninety angular blocks of 2-6 m radius plus six rounded coral heads,
+		# on rippled sand.
+		SiteTerrainBaker.BoulderField.new(BOULDER_FIELD_CENTRE, BOULDER_FIELD_RADIUS_M, 301, 90, 2.0, 6.0, 0.7, 0.85),
+		SiteTerrainBaker.BoulderField.new(Vector2(-100.0, 130.0), 40.0, 302, 6, 6.0, 9.0, 0.6, 0.0),
 		SiteTerrainBaker.Ripples.new(Vector2(1.0, 0.3), 3.5, 0.3, BOULDER_FIELD_CENTRE, 110.0),
 
 		# SE — les Éperons: two pitted rock spurs 7-8 m high framing a 20 m sand channel
 		# littered with small blocks.
 		SiteTerrainBaker.Outcrop.new(WEST_SPUR_CENTRE, 30.0, 60.0, 8.0, 5.0, 401, 7.0, true),
 		SiteTerrainBaker.Outcrop.new(EAST_SPUR_CENTRE, 30.0, 55.0, 7.0, 5.0, 402, 7.0, true),
-		SiteTerrainBaker.BoulderField.new(CHANNEL_CENTRE, 50.0, 403, 12, 1.5, 3.0, 0.8),
+		SiteTerrainBaker.BoulderField.new(CHANNEL_CENTRE, 50.0, 403, 12, 1.5, 3.0, 0.8, 1.0),
 
 		# S — la Corniche: sinuous 22 m escarpment, then a 42 degree slope to the 80 m floor,
 		# with rounded knolls on the slope.
@@ -97,8 +98,13 @@ static func recipe() -> SiteTerrainBaker.Recipe:
 	result.modifiers = [
 		# Crumbly rock faces, smooth sand. Before the strata so the ledges stay flat.
 		SiteTerrainBaker.SlopeRoughness.new(2028, 0.12, 0.7, 20.0, 40.0),
-		# Strata on the north-west quarter only (mesa flanks and canyon walls): 4 m ledges.
-		SiteTerrainBaker.Strata.new(4.0, 0.85, 30.0, 50.0, Rect2(-256.0, -256.0, 256.0, 256.0)),
+		# Strata on the north-west quarter only (mesa flanks and canyon walls): beds of
+		# about 4 m, dipping 1.3 degrees to the south-east, bent by 2.4 m of slow warp,
+		# each bed with its own thickness and hardness.
+		SiteTerrainBaker.Strata.new(
+			4.0, 0.85, 30.0, 50.0, Rect2(-256.0, -256.0, 256.0, 256.0),
+			Vector2(0.015, 0.018), 105, 2.4, 0.5
+		),
 	]
 	return result
 

tools/sites/zone1_vigilant_recipe.gd

@@ -75,7 +75,10 @@ static func recipe() -> SiteTerrainBaker.Recipe:
 		# Crumbly rock on every slope steeper than the sand rule, smooth sand elsewhere.
 		# Before the strata so the ledges stay flat.
 		SiteTerrainBaker.SlopeRoughness.new(1935, 0.12, 0.7, 20.0, 40.0),
-		SiteTerrainBaker.Strata.new(4.0, 0.85, 25.0, 40.0, Rect2(-384.0, STRATA_START_Z, 768.0, 384.0 - STRATA_START_Z)),
+		SiteTerrainBaker.Strata.new(
+			4.0, 0.85, 25.0, 40.0, Rect2(-384.0, STRATA_START_Z, 768.0, 384.0 - STRATA_START_Z),
+			Vector2(0.008, 0.0), 1938, 2.0, 0.45
+		),
 	]
 	result.material_patches = [
 		SiteTerrainBaker.MaterialPatch.new(Vector2(140.0, 40.0), 28.0, SonarMaterial.Value.ROCK),