NoRaincheck · January 27, 2026 00:29
diff --git a/diamond_square.py b/diamond_square.py
 """
 Diamond-Square Terrain Generator

 A Python implementation of the Diamond-Square algorithm for procedural fractal
 terrain generation, with extensions for constrained generation, seamless tile
 stitching, and water ratio control.

 THE DIAMOND-SQUARE ALGORITHM
 ============================

 The Diamond-Square algorithm generates realistic fractal terrain through an
 iterative subdivision process:

 1. Initialize corners: Set the four corners of a square grid with initial values
 2. Diamond step: For each square, set the center point to the average of the
   four corners plus a random displacement
 3. Square step: For each diamond, set the center point to the average of the
   four adjacent points plus a random displacement
 4. Repeat: Halve the step size, reduce the random displacement by the roughness
   factor, and repeat until all points are filled

 The grid size must be 2^n + 1 (e.g., 129, 257, 513) to ensure proper subdivision.

 Roughness Parameter:
 - Low roughness (0.2-0.4): Smooth, rolling hills
 - Medium roughness (0.5): Balanced terrain (default)
 - High roughness (0.7-0.9): Jagged, mountainous terrain

 EXTENSIONS
 ==========

 Constrained Generation:
 - Fixed corners: Set specific corner values via the `corners` parameter
 - Arbitrary constraints: Pre-set any (y, x) position via the `constraints` dict
 - Constrained points are preserved during generation

 Border Seeding:
 - The `seed_border` option pre-seeds all border edges with low values
 - Creates island-like terrain where edges naturally become water

 Seamless Tile Stitching:
 - The TerrainGrid class enables infinite terrain generation
 - Tiles share edge values with their neighbors
 - New tiles are generated with constraints from already-generated adjacent tiles
 - Supports deterministic seeding per tile for reproducibility

 Grid Modes (--grid vs --grid2):
 - --grid: "True to arguments" mode. Generates exactly what you specify with no
  automatic adjustments. The output dimensions and water ratio are directly
  determined by your parameters without any correction.
  
 - --grid2: "Best effort/smart" mode. Adds automatic postprocessing to improve results:
  * Corner smoothing: Blends the top-left corner toward water with noisy falloff,
    reducing visible tile seams at the origin.
  * When combined with --border: Enforces water borders by iteratively increasing
    water ratio until all edges are water, then applies smart-crop to recenter
    land and minimize water margins.
  * Terrain stacking: If water ratio diverges too far from target (>10%), blends
    multiple terrain generations using optimal EMA weights to recover the target.
  
  Use --grid for predictable output; use --grid2 when you want the generator to
  "do its best" to produce island-like terrain with water borders and reasonable
  water ratios.

 WATER RATIO CONTROL
 ===================

 Two methods are available to achieve a target water ratio:

 Piecewise Method (Exact):
 - Finds the value at the target percentile
 - Remaps values below to [0, WATER_THRESHOLD]
 - Remaps values above to [WATER_THRESHOLD, 1.0]
 - Guarantees exactly the requested water ratio

 Gamma Method (Default, Approximate):
 - Uses power transform (terrain^γ) with binary search
 - Higher gamma pushes values toward 0 (more water)
 - Lower gamma pushes values toward 1 (less water)
 - Produces smoother results but may not hit exact target (typically within 0.1%)
 """

 import argparse
 import numpy as np
 from PIL import Image
 from scipy.ndimage import zoom

 # Water threshold: values below this are considered water (based on colormap)
 # Coast water ends at 50/255 in the colormap
 WATER_THRESHOLD = 50 / 255  # ~0.196

 # Border enforcement parameters
 BORDER_ENFORCEMENT_MAX_RATIO = 0.95
 BORDER_ENFORCEMENT_INCREMENT = 0.05

 # Terrain stacking parameters (for --grid2 --border mode)
 STACKING_MAX_ITERATIONS = 5
 STACKING_WATER_RATIO_TOLERANCE = 0.10  # 10%

 # Corner indices for consistent ordering: (y, x) positions
 CORNER_POSITIONS = [
    (0, 0),      # top-left
    (0, -1),     # top-right (use size-1 when indexing)
    (-1, 0),     # bottom-left
    (-1, -1),    # bottom-right
 ]


 def diamond_square(
    n: int,
    roughness: float = 0.5,
    seed: int | None = None,
    corners: tuple[float, float, float, float] | None = None,
    constraints: dict[tuple[int, int], float] | None = None,
    normalize: bool = True,
    seed_border: bool | float = False,
 ) -> np.ndarray:
    """
    Generate terrain using the Diamond-Square algorithm.
    
    Args:
        n: Power of 2 for grid size (grid will be 2^n + 1 on each side)
        roughness: Decay factor for random displacement (0.0-1.0)
        seed: Optional random seed for reproducibility
        corners: Optional fixed corner values (top-left, top-right, bottom-left, bottom-right)
        constraints: Dictionary of (y, x) -> value for pre-set values that won't be modified
        normalize: Whether to normalize output to 0-1 range (default True)
        seed_border: Pre-seed all border edges with a low value.
                     If True: use -2*roughness. If float: use that value. If False: no seeding.
    
    Returns:
        2D numpy array of terrain heights (normalized to 0-1 if normalize=True)
    """
    rng = np.random.RandomState(seed)
    size = (1 << n) + 1  # 2^n + 1
    terrain = np.full((size, size), np.nan)
    
    # Apply constraints if provided
    if constraints:
        for (y, x), value in constraints.items():
            if 0 <= y < size and 0 <= x < size:
                terrain[y, x] = value
    
    # Apply border seeding if requested (only to positions not already set)
    if seed_border is not False:
        border_val = -2 * roughness if seed_border is True else seed_border
        _seed_terrain_border(terrain, border_val)
    
    # Initialize corners
    corner_coords = [(0, 0), (0, size - 1), (size - 1, 0), (size - 1, size - 1)]
    corner_values = corners if corners else tuple(rng.normal(scale=roughness) for _ in range(4))
    
    for (y, x), value in zip(corner_coords, corner_values):
        if np.isnan(terrain[y, x]):
            terrain[y, x] = value
    
    # Diamond-square iteration
    step_size = size - 1
    scale = 1.0
    
    while step_size > 1:
        half = step_size // 2
        _diamond_step(terrain, size, step_size, half, scale, rng)
        _square_step(terrain, size, step_size, half, scale, rng)
        step_size = half
        scale *= roughness
    
    if normalize:
        terrain = _normalize_terrain(terrain)
    
    return terrain


 def _seed_terrain_border(terrain: np.ndarray, border_val: float) -> None:
    """Seed all border edges with a value (in-place)."""
    size = terrain.shape[0]
    for i in range(size):
        for y, x in [(0, i), (size - 1, i), (i, 0), (i, size - 1)]:
            if np.isnan(terrain[y, x]):
                terrain[y, x] = border_val


 def _diamond_step(
    terrain: np.ndarray,
    size: int,
    step_size: int,
    half: int,
    scale: float,
    rng: np.random.RandomState,
 ) -> None:
    """Compute center of each square (diamond step)."""
    for y in range(half, size, step_size):
        for x in range(half, size, step_size):
            if np.isnan(terrain[y, x]):
                avg = (
                    terrain[y - half, x - half] +
                    terrain[y - half, x + half] +
                    terrain[y + half, x - half] +
                    terrain[y + half, x + half]
                ) / 4.0
                terrain[y, x] = avg + rng.normal(scale=scale)


 def _square_step(
    terrain: np.ndarray,
    size: int,
    step_size: int,
    half: int,
    scale: float,
    rng: np.random.RandomState,
 ) -> None:
    """Compute center of each diamond (square step)."""
    for y in range(0, size, half):
        for x in range((y + half) % step_size, size, step_size):
            if np.isnan(terrain[y, x]):
                neighbors = [
                    (y - half, x), (y + half, x),
                    (y, x - half), (y, x + half),
                ]
                total = sum(
                    terrain[ny, nx]
                    for ny, nx in neighbors
                    if 0 <= ny < size and 0 <= nx < size
                )
                count = sum(
                    1 for ny, nx in neighbors
                    if 0 <= ny < size and 0 <= nx < size
                )
                terrain[y, x] = (total / count) + rng.normal(scale=scale)


 def _normalize_terrain(terrain: np.ndarray) -> np.ndarray:
    """Normalize terrain values to 0-1 range."""
    min_val = np.nanmin(terrain)
    max_val = np.nanmax(terrain)
    if max_val > min_val:
        return (terrain - min_val) / (max_val - min_val)
    return np.zeros_like(terrain)


 def rescale_for_water_ratio(
    terrain: np.ndarray,
    target_ratio: float,
    method: str = "gamma",
 ) -> np.ndarray:
    """
    Rescale terrain to achieve a target water ratio.
    
    Args:
        terrain: 2D array of terrain heights (normalized to 0-1)
        target_ratio: Target proportion of water (0.0-1.0)
        method: Rescaling method - "gamma" or "piecewise"
    
    Returns:
        Rescaled terrain with approximately target_ratio values below WATER_THRESHOLD
    """
    rescalers = {
        "gamma": _rescale_gamma,
        "piecewise": _rescale_piecewise,
    }
    if method not in rescalers:
        raise ValueError(f"Unknown rescaling method: {method}")
    return rescalers[method](terrain, target_ratio)


 def _rescale_gamma(terrain: np.ndarray, target_ratio: float, max_iterations: int = 50) -> np.ndarray:
    """
    Gamma/power transform rescaling.
    
    Binary search for gamma such that terrain^gamma achieves the target water ratio.
    Higher gamma pushes values toward 0 (more water), lower gamma toward 1 (less water).
    """
    gamma_low, gamma_high = 0.1, 10.0
    gamma = (gamma_low + gamma_high) / 2
    
    for _ in range(max_iterations):
        gamma = (gamma_low + gamma_high) / 2
        transformed = np.power(terrain, gamma)
        current_ratio = np.mean(transformed < WATER_THRESHOLD)
        
        if abs(current_ratio - target_ratio) < 0.001:
            break
        
        if current_ratio < target_ratio:
            gamma_low = gamma  # Need more water -> higher gamma
        else:
            gamma_high = gamma  # Need less water -> lower gamma
    
    return np.power(terrain, gamma)


 def _rescale_piecewise(terrain: np.ndarray, target_ratio: float) -> np.ndarray:
    """
    Piecewise linear rescaling.
    
    Find the value at target percentile, then rescale:
    - Below percentile: stretch/compress to fill [0, WATER_THRESHOLD]
    - Above percentile: stretch/compress to fill [WATER_THRESHOLD, 1]
    """
    percentile_value = np.percentile(terrain, target_ratio * 100)
    min_val, max_val = terrain.min(), terrain.max()
    
    if max_val <= min_val or percentile_value <= min_val or percentile_value >= max_val:
        return terrain.copy()
    
    result = np.empty_like(terrain)
    below_mask = terrain <= percentile_value
    
    # Map [min_val, percentile_value] -> [0, WATER_THRESHOLD]
    below_range = percentile_value - min_val
    result[below_mask] = (
        (terrain[below_mask] - min_val) / below_range * WATER_THRESHOLD
        if below_range > 0 else 0
    )
    
    # Map [percentile_value, max_val] -> [WATER_THRESHOLD, 1]
    above_range = max_val - percentile_value
    above_mask = ~below_mask
    result[above_mask] = (
        WATER_THRESHOLD + (terrain[above_mask] - percentile_value) / above_range * (1 - WATER_THRESHOLD)
        if above_range > 0 else 1.0
    )
    
    return np.clip(result, 0, 1)


 def blur_water_regions(terrain: np.ndarray) -> np.ndarray:
    """
    Apply gaussian blur to water regions only.
    
    Args:
        terrain: 2D array of terrain heights (normalized to 0-1)
    
    Returns:
        Terrain with blurred water regions (smooths coastlines)
    """
    from scipy.ndimage import gaussian_filter
    
    # Auto-scale sigma based on image size (0.5% of smaller dimension)
    sigma = max(1, min(terrain.shape) * 0.005)
    
    water_mask = terrain < WATER_THRESHOLD
    blurred = gaussian_filter(terrain, sigma=sigma)
    
    result = terrain.copy()
    result[water_mask] = blurred[water_mask]
    return result


 def smart_crop(terrain: np.ndarray, margin_ratio: float = 0.05) -> np.ndarray:
    """
    Intelligently crop terrain to minimize water border while keeping land centered.
    
    Args:
        terrain: 2D array of terrain heights (normalized to 0-1)
        margin_ratio: Margin to add around land bounding box (as ratio of bbox size)
    
    Returns:
        Cropped and resized terrain at original dimensions
    """
    original_shape = terrain.shape
    land_mask = terrain >= WATER_THRESHOLD
    
    # Find land bounding box
    rows = np.any(land_mask, axis=1)
    cols = np.any(land_mask, axis=0)
    if not rows.any() or not cols.any():
        return terrain  # No land found
    
    row_min, row_max = np.where(rows)[0][[0, -1]]
    col_min, col_max = np.where(cols)[0][[0, -1]]
    
    # Add margin
    height, width = row_max - row_min, col_max - col_min
    margin_y, margin_x = int(height * margin_ratio), int(width * margin_ratio)
    
    row_min = max(0, row_min - margin_y)
    row_max = min(terrain.shape[0] - 1, row_max + margin_y)
    col_min = max(0, col_min - margin_x)
    col_max = min(terrain.shape[1] - 1, col_max + margin_x)
    
    # Crop and resize back to original dimensions
    cropped = terrain[row_min:row_max + 1, col_min:col_max + 1]
    zoom_factors = (original_shape[0] / cropped.shape[0], original_shape[1] / cropped.shape[1])
    return zoom(cropped, zoom_factors, order=1)


 def has_water_border(terrain: np.ndarray, threshold: float = 1.0) -> bool:
    """
    Check if all four edges of the terrain are entirely water.
    
    Args:
        terrain: 2D array of terrain heights (normalized to 0-1)
        threshold: Minimum fraction of border pixels that must be water (default 1.0 = all)
    
    Returns:
        True if all border pixels are water (value < WATER_THRESHOLD)
    """
    top_edge = terrain[0, :]
    bottom_edge = terrain[-1, :]
    left_edge = terrain[1:-1, 0]
    right_edge = terrain[1:-1, -1]
    
    all_border = np.concatenate([top_edge, bottom_edge, left_edge, right_edge])
    return np.mean(all_border < WATER_THRESHOLD) >= threshold


 def enforce_water_border(
    terrain: np.ndarray,
    water_method: str,
    initial_ratio: float | None = None,
 ) -> tuple[np.ndarray, bool, float]:
    """
    Iteratively increase water ratio until all borders are water.
    
    Args:
        terrain: 2D array of terrain heights (normalized to 0-1)
        water_method: Method for water ratio rescaling
        initial_ratio: Starting water ratio (defaults to current ratio)
    
    Returns:
        Tuple of (modified terrain, success, final_ratio)
    """
    if has_water_border(terrain):
        return terrain, True, np.mean(terrain < WATER_THRESHOLD)
    
    current_ratio = initial_ratio if initial_ratio is not None else np.mean(terrain < WATER_THRESHOLD)
    current_ratio = max(current_ratio, 0.1)
    
    while not has_water_border(terrain) and current_ratio < BORDER_ENFORCEMENT_MAX_RATIO:
        current_ratio = min(current_ratio + BORDER_ENFORCEMENT_INCREMENT, BORDER_ENFORCEMENT_MAX_RATIO)
        terrain = rescale_for_water_ratio(terrain, current_ratio, water_method)
    
    return terrain, has_water_border(terrain), current_ratio


 def stack_terrains_for_water_ratio(
    terrain: np.ndarray,
    target_ratio: float,
    args: "argparse.Namespace",
    grid_w: int,
    grid_h: int,
 ) -> tuple[np.ndarray, int]:
    """
    Iteratively stack terrains to bring water ratio closer to target.
    
    When --grid2 --border causes water ratio to exceed target significantly,
    generate additional terrains and combine them using weighted averaging
    until the water ratio is within tolerance.
    
    The stacking approach works by averaging multiple terrains, which tends to
    reduce extreme values. Since borders need to be water, we track the best
    result that still maintains water borders and is closest to target.
    
    Args:
        terrain: Initial terrain (already processed with enforce_water_border)
        target_ratio: Target water ratio (0.0-1.0)
        args: Command line arguments for terrain generation
        grid_w: Grid width in tiles
        grid_h: Grid height in tiles
    
    Returns:
        Tuple of (final terrain, number of stacking iterations performed)
    """
    current_ratio = np.mean(terrain < WATER_THRESHOLD)
    
    # Check if already within tolerance
    if abs(current_ratio - target_ratio) <= STACKING_WATER_RATIO_TOLERANCE:
        return terrain, 0
    
    # Only stack if water ratio is too high (common case with border enforcement)
    if current_ratio <= target_ratio:
        return terrain, 0
    
    # Track the best terrain that still has water borders
    best_terrain = terrain.copy()
    best_ratio = current_ratio
    best_error = abs(current_ratio - target_ratio)
    
    # Running stacked terrain for EMA blending
    stacked = terrain.copy()
    
    base_seed = args.seed if args.seed is not None else 0
    iterations = 0
    
    for iteration in range(1, STACKING_MAX_ITERATIONS + 1):
        # Generate new terrain with different seed
        new_seed = base_seed + iteration * 1000
        
        grid = TerrainGrid(
            args.power,
            args.roughness,
            new_seed,
            seed_global_border=args.border,
            smooth_corner=args.grid2 is not None,
            target_water_ratio=target_ratio,
            water_method=args.water_method,
        )
        new_terrain = grid.get_region(0, 0, grid_w, grid_h)
        
        # Enforce water border on the new terrain before blending
        new_terrain, _, _ = enforce_water_border(new_terrain, args.water_method, target_ratio)
        
        # Use scipy.optimize to find the optimal alpha for blending
        from scipy.optimize import minimize_scalar
        
        def objective(alpha: float) -> float:
            """Compute water ratio error for a given blend alpha."""
            blended = (1 - alpha) * stacked + alpha * new_terrain
            blended = _normalize_terrain(blended)
            blended = rescale_for_water_ratio(blended, target_ratio, args.water_method)
            blended, borders_ok, _ = enforce_water_border(blended, args.water_method, target_ratio)
            ratio = np.mean(blended < WATER_THRESHOLD)
            # Penalize if borders are broken
            penalty = 0.0 if borders_ok else 1.0
            return abs(ratio - target_ratio) + penalty
        
        # Find optimal alpha in (0.01, 0.99) - avoid extremes
        result = minimize_scalar(objective, bounds=(0.01, 0.99), method='bounded')
        optimal_alpha = result.x
        
        # Apply the optimal blend
        stacked = (1 - optimal_alpha) * stacked + optimal_alpha * new_terrain
        
        # Re-normalize after combining
        stacked = _normalize_terrain(stacked)
        
        # Re-apply water ratio rescaling toward target
        stacked = rescale_for_water_ratio(stacked, target_ratio, args.water_method)
        
        # Re-enforce water borders (normalization/rescaling may have broken them)
        stacked, borders_ok, enforced_ratio = enforce_water_border(stacked, args.water_method, target_ratio)
        
        iterations = iteration
        new_ratio = np.mean(stacked < WATER_THRESHOLD)
        new_error = abs(new_ratio - target_ratio)
        
        update_message = ""
        # If this result is better AND maintains water borders, save it
        if borders_ok and new_error < best_error:
            update_message = "...Updating Terrain"
            best_terrain = stacked.copy()
            best_ratio = new_ratio
            best_error = new_error
        
        print(f"  Stacking iteration {iteration}: water ratio {current_ratio:.3f} -> {new_ratio:.3f} (alpha={optimal_alpha:.3f}){update_message}")
        stacked = best_terrain.copy()
        current_ratio = np.mean(stacked < WATER_THRESHOLD)
        current_error = abs(current_ratio - target_ratio)
        
        # Check if we've reached tolerance (with water borders intact)
        if borders_ok and new_error <= STACKING_WATER_RATIO_TOLERANCE:
            return stacked, iterations
    
    # Return the best result we found
    print(f"  Best result: water ratio {best_ratio:.3f} (error: {best_error:.3f})")
    return best_terrain, iterations


 def remove_small_islands(terrain: np.ndarray, size_threshold: float = 0.05) -> tuple[np.ndarray, int]:
    """
    Remove small land masses by converting them to water.
    
    Uses connected component labeling to identify separate land regions,
    then removes any that are smaller than a threshold percentage of the
    largest landmass.
    
    Args:
        terrain: 2D array of terrain heights (normalized to 0-1)
        size_threshold: Islands smaller than this fraction of the largest landmass are removed
    
    Returns:
        (modified terrain, number of islands removed)
    """
    from scipy.ndimage import label
    
    land_mask = terrain >= WATER_THRESHOLD
    labeled_array, num_features = label(land_mask)
    
    if num_features <= 1:
        return terrain, 0
    
    # Find size of each region
    region_sizes = {
        region_id: np.sum(labeled_array == region_id)
        for region_id in range(1, num_features + 1)
    }
    
    min_size = max(region_sizes.values()) * size_threshold
    islands_to_remove = [rid for rid, size in region_sizes.items() if size < min_size]
    
    if not islands_to_remove:
        return terrain, 0
    
    result = terrain.copy()
    for region_id in islands_to_remove:
        result[labeled_array == region_id] = WATER_THRESHOLD * 0.5
    
    return result, len(islands_to_remove)


 class TerrainGrid:
    """
    Manages infinite terrain generation via constrained tile stitching.
    
    Uses a sliding window approach where new tiles are generated with edge
    constraints from already-generated neighboring tiles.
    """
    
    # Neighbor offsets: (dx, dy) -> (edge_func for neighbor, edge_func for self)
    # edge_func takes (tile, size) and returns list of ((y, x), value) constraints
    NEIGHBOR_EDGES = {
        (-1, 0): (  # left neighbor
            lambda t, s: [(i, s - 1) for i in range(s)],  # neighbor's right edge
            lambda s: [(i, 0) for i in range(s)],         # our left edge
        ),
        (1, 0): (   # right neighbor
            lambda t, s: [(i, 0) for i in range(s)],      # neighbor's left edge
            lambda s: [(i, s - 1) for i in range(s)],     # our right edge
        ),
        (0, -1): (  # top neighbor
            lambda t, s: [(s - 1, i) for i in range(s)],  # neighbor's bottom edge
            lambda s: [(0, i) for i in range(s)],         # our top edge
        ),
        (0, 1): (   # bottom neighbor
            lambda t, s: [(0, i) for i in range(s)],      # neighbor's top edge
            lambda s: [(s - 1, i) for i in range(s)],     # our bottom edge
        ),
    }
    
    DIAGONAL_CORNERS = {
        (-1, -1): ((lambda s: s - 1, lambda s: s - 1), (0, 0)),           # top-left
        (1, -1):  ((lambda s: s - 1, lambda s: 0), (0, lambda s: s - 1)), # top-right
        (-1, 1):  ((lambda s: 0, lambda s: s - 1), (lambda s: s - 1, 0)), # bottom-left
        (1, 1):   ((lambda s: 0, lambda s: 0), (lambda s: s - 1, lambda s: s - 1)), # bottom-right
    }
    
    def __init__(
        self,
        tile_n: int,
        roughness: float = 0.5,
        seed: int | None = None,
        seed_global_border: bool = False,
        smooth_corner: bool = False,
        target_water_ratio: float | None = None,
        water_method: str = "percentile",
        polyak_alpha: float = 0.1,
    ):
        """
        Initialize the terrain grid.
        
        Args:
            tile_n: Power of 2 for tile size (each tile will be 2^n + 1 on each side)
            roughness: Decay factor for random displacement (0.0-1.0)
            seed: Base random seed for reproducibility
            seed_global_border: If True, seed the outer edges of the region with low values
            smooth_corner: If True, postprocess to smooth top-left corner toward running min
            target_water_ratio: Target proportion of water (0.0-1.0), or None to disable
            water_method: Method for water ratio rescaling ("percentile", "gamma", "piecewise")
            polyak_alpha: Exponential moving average factor for running statistics (0.0-1.0)
        """
        self.tile_n = tile_n
        self.tile_size = (1 << tile_n) + 1
        self.roughness = roughness
        self.base_seed = seed if seed is not None else np.random.RandomState().randint(0, 2**31)
        self.tiles: dict[tuple[int, int], np.ndarray] = {}
        self.seed_global_border = seed_global_border
        self.smooth_corner = smooth_corner
        self._running_min: float = -2 * roughness
        
        # Water ratio settings
        self._target_water_ratio = target_water_ratio
        self._water_method = water_method
        self._polyak_alpha = polyak_alpha
        
        # Polyak averaging for distribution estimation
        self._percentile_points = np.arange(5, 100, 5)
        self._running_percentiles: np.ndarray | None = None
        self._running_min_raw: float | None = None
        self._running_max_raw: float | None = None
        self._tile_count: int = 0
    
    def _get_tile_seed(self, x: int, y: int) -> int:
        """Generate a deterministic seed for a tile position."""
        return hash((self.base_seed, x, y)) & 0x7FFFFFFF
    
    def _update_running_stats(self, tile: np.ndarray) -> None:
        """Update running statistics using Polyak (exponential moving) averaging."""
        tile_percentiles = np.percentile(tile.flatten(), self._percentile_points)
        tile_min, tile_max = float(tile.min()), float(tile.max())
        
        if self._running_percentiles is None:
            self._running_percentiles = tile_percentiles
            self._running_min_raw = tile_min
            self._running_max_raw = tile_max
        else:
            alpha = self._polyak_alpha
            self._running_percentiles = (1 - alpha) * self._running_percentiles + alpha * tile_percentiles
            self._running_min_raw = (1 - alpha) * self._running_min_raw + alpha * tile_min
            self._running_max_raw = (1 - alpha) * self._running_max_raw + alpha * tile_max
        
        self._tile_count += 1
    
    def _get_neighbor_constraints(self, x: int, y: int) -> dict[tuple[int, int], float]:
        """Build constraints from neighboring tiles."""
        constraints: dict[tuple[int, int], float] = {}
        size = self.tile_size
        
        # Edge neighbors
        for (dx, dy), (neighbor_edge, self_edge) in self.NEIGHBOR_EDGES.items():
            neighbor_pos = (x + dx, y + dy)
            if neighbor_pos in self.tiles:
                neighbor_tile = self.tiles[neighbor_pos]
                for (ny, nx), (sy, sx) in zip(neighbor_edge(None, size), self_edge(size)):
                    constraints[(sy, sx)] = neighbor_tile[ny, nx]
        
        # Diagonal corners
        for (dx, dy), (neighbor_idx, self_idx) in self.DIAGONAL_CORNERS.items():
            neighbor_pos = (x + dx, y + dy)
            if neighbor_pos in self.tiles:
                ny = neighbor_idx[0](size) if callable(neighbor_idx[0]) else neighbor_idx[0]
                nx = neighbor_idx[1](size) if callable(neighbor_idx[1]) else neighbor_idx[1]
                sy = self_idx[0](size) if callable(self_idx[0]) else self_idx[0]
                sx = self_idx[1](size) if callable(self_idx[1]) else self_idx[1]
                constraints[(sy, sx)] = self.tiles[neighbor_pos][ny, nx]
        
        return constraints
    
    def _apply_border_constraints(
        self,
        constraints: dict[tuple[int, int], float],
        x: int,
        y: int,
        global_bounds: tuple[int, int, int, int],
    ) -> None:
        """Apply global border constraints for tiles at region boundaries."""
        min_x, min_y, max_x, max_y = global_bounds
        size = self.tile_size
        border_val = self._running_min
        
        border_edges = [
            (x == min_x, [(i, 0) for i in range(size)]),         # left edge
            (x == max_x, [(i, size - 1) for i in range(size)]),  # right edge
            (y == min_y, [(0, i) for i in range(size)]),         # top edge
            (y == max_y, [(size - 1, i) for i in range(size)]),  # bottom edge
        ]
        
        for is_boundary, positions in border_edges:
            if is_boundary:
                for pos in positions:
                    if pos not in constraints:
                        constraints[pos] = border_val
    
    def get_tile(
        self,
        x: int,
        y: int,
        global_bounds: tuple[int, int, int, int] | None = None,
    ) -> np.ndarray:
        """
        Get or generate a tile at the given grid position.
        
        Args:
            x: Tile x-coordinate in the grid
            y: Tile y-coordinate in the grid
            global_bounds: Optional (min_x, min_y, max_x, max_y) defining the region bounds.
        
        Returns:
            2D numpy array of terrain heights for this tile
        """
        if (x, y) in self.tiles:
            return self.tiles[(x, y)]
        
        constraints = self._get_neighbor_constraints(x, y)
        
        if self.seed_global_border and global_bounds is not None:
            self._apply_border_constraints(constraints, x, y, global_bounds)
        
        tile = diamond_square(
            n=self.tile_n,
            roughness=self.roughness,
            seed=self._get_tile_seed(x, y),
            constraints=constraints if constraints else None,
            normalize=False,
        )
        
        # Update running statistics
        tile_min = np.min(tile)
        if tile_min < self._running_min:
            self._running_min = tile_min
        self._update_running_stats(tile)
        
        self.tiles[(x, y)] = tile
        return tile
    
    def _apply_corner_smoothing(self, result: np.ndarray) -> None:
        """Apply corner smoothing postprocessing (for --grid2 mode)."""
        size = self.tile_size
        blend_into_neighbors = size // 4
        affected_rows = size + blend_into_neighbors
        affected_cols = size + blend_into_neighbors
        max_radius = np.sqrt(affected_rows**2 + affected_cols**2)
        target_min = self._running_min
        
        rng = np.random.RandomState(self.base_seed + 999)
        
        for row in range(min(affected_rows, result.shape[0])):
            for col in range(min(affected_cols, result.shape[1])):
                dist = np.sqrt(row**2 + col**2)
                falloff = max(0, 1 - dist / max_radius)
                
                if falloff > 0:
                    noise = rng.normal(0, self.roughness * 0.5)
                    noisy_target = target_min + abs(noise)
                    blend_strength = falloff ** 2
                    blended = result[row, col] * (1 - blend_strength) + noisy_target * blend_strength
                    result[row, col] = (result[row, col] * 2 + blended) / 3
    
    def get_region(self, grid_x: int, grid_y: int, width: int, height: int) -> np.ndarray:
        """
        Get a stitched region spanning multiple tiles.
        
        Args:
            grid_x: Starting tile x-coordinate
            grid_y: Starting tile y-coordinate
            width: Number of tiles wide
            height: Number of tiles tall
        
        Returns:
            2D numpy array of the combined terrain
        """
        tile_inner = self.tile_size - 1
        output_height = height * tile_inner + 1
        output_width = width * tile_inner + 1
        result = np.zeros((output_height, output_width))
        
        global_bounds = (
            (grid_x, grid_y, grid_x + width - 1, grid_y + height - 1)
            if self.seed_global_border else None
        )
        
        # Generate tiles in order to ensure constraints are available
        for ty in range(height):
            for tx in range(width):
                tile = self.get_tile(grid_x + tx, grid_y + ty, global_bounds=global_bounds)
                out_y, out_x = ty * tile_inner, tx * tile_inner
                result[out_y:out_y + self.tile_size, out_x:out_x + self.tile_size] = tile
        
        if self.smooth_corner:
            self._apply_corner_smoothing(result)
        
        result = _normalize_terrain(result)
        
        if self._target_water_ratio is not None:
            result = rescale_for_water_ratio(result, self._target_water_ratio, self._water_method)
        
        return result
    
    def clear_cache(self) -> None:
        """Clear the tile cache and running statistics."""
        self.tiles.clear()
        self._running_percentiles = None
        self._running_min_raw = None
        self._running_max_raw = None
        self._tile_count = 0


 def create_terrain_colormap() -> np.ndarray:
    """
    Create a terrain colormap with natural colors.
    
    Returns:
        256x3 numpy array of RGB values
    """
    colormap = np.zeros((256, 3), dtype=np.uint8)
    
    stops = [
        (0,    (28,  60,  112)),   # Deep water
        (35,   (47,  102, 144)),   # Shallow water
        (50,   (70,  130, 160)),   # Coast water
        (55,   (194, 178, 128)),   # Beach
        (65,   (85,  128, 60)),    # Lowland
        (100,  (64,  105, 48)),    # Forest
        (140,  (95,  95,  75)),    # Highland
        (180,  (128, 118, 105)),   # Mountain
        (210,  (168, 160, 150)),   # High mountain
        (235,  (210, 210, 215)),   # Alpine
        (255,  (250, 250, 255)),   # Snow peak
    ]
    
    for i, (start_idx, start_color) in enumerate(stops[:-1]):
        end_idx, end_color = stops[i + 1]
        for j in range(start_idx, end_idx + 1):
            t = (j - start_idx) / (end_idx - start_idx) if end_idx > start_idx else 0
            colormap[j] = [
                int(start_color[c] + t * (end_color[c] - start_color[c]))
                for c in range(3)
            ]
    
    return colormap


 def save_terrain_image(terrain: np.ndarray, filename: str) -> None:
    """
    Save terrain data as a colored terrain PNG image.
    
    Args:
        terrain: 2D array of terrain heights (normalized to 0-1)
        filename: Output filename
    """
    terrain_uint8 = (terrain * 255).astype(np.uint8)
    colormap = create_terrain_colormap()
    colored = colormap[terrain_uint8]
    
    Image.fromarray(colored, mode='RGB').save(filename)
    print(f"Saved terrain image to: {filename}")


 def parse_grid_spec(spec: str) -> tuple[int, int] | None:
    """Parse a grid specification string like '3x3' into (width, height)."""
    try:
        return tuple(map(int, spec.lower().split('x')))
    except ValueError:
        return None


 def generate_single_terrain(args: argparse.Namespace) -> np.ndarray:
    """Generate a single-tile terrain based on command line arguments."""
    terrain = diamond_square(args.power, args.roughness, args.seed, seed_border=args.border)
    if args.water_ratio is not None:
        terrain = rescale_for_water_ratio(terrain, args.water_ratio, args.water_method)
    return terrain


 def generate_grid_terrain(args: argparse.Namespace, grid_w: int, grid_h: int) -> np.ndarray:
    """Generate a multi-tile grid terrain based on command line arguments."""
    grid = TerrainGrid(
        args.power,
        args.roughness,
        args.seed,
        seed_global_border=args.border,
        smooth_corner=args.grid2 is not None,
        target_water_ratio=args.water_ratio,
        water_method=args.water_method,
    )
    return grid.get_region(0, 0, grid_w, grid_h)


 def apply_postprocessing(terrain: np.ndarray, args: argparse.Namespace) -> np.ndarray:
    """Apply postprocessing steps (smart crop, water blur, island removal)."""
    if not args.smart_crop:
        if args.water_blur:
            terrain = blur_water_regions(terrain)
            print("Applied gaussian blur to water regions")
        return terrain
    
    pre_crop_ratio = np.mean(terrain < WATER_THRESHOLD)
    
    # Remove small islands if water ratio diverged significantly
    if args.water_ratio is not None and abs(pre_crop_ratio - args.water_ratio) >= 0.25:
        terrain, removed = remove_small_islands(terrain)
        if removed > 0:
            new_ratio = np.mean(terrain < WATER_THRESHOLD)
            print(f"Removed {removed} small island(s), water ratio: {pre_crop_ratio:.3f} -> {new_ratio:.3f}")
            pre_crop_ratio = new_ratio
    
    terrain = smart_crop(terrain)
    post_crop_ratio = np.mean(terrain < WATER_THRESHOLD)
    delta = post_crop_ratio - pre_crop_ratio
    print(f"Applied smart crop to minimize water border")
    print(f"Water ratio change: {pre_crop_ratio:.3f} -> {post_crop_ratio:.3f} ({delta:+.3f})")
    
    terrain = blur_water_regions(terrain)
    print("Applied gaussian blur to water regions")
    
    return terrain


 def main():
    parser = argparse.ArgumentParser(
        description="Generate fractal terrain using the Diamond-Square algorithm"
    )
    parser.add_argument(
        "-n", "--power", type=int, default=9,
        help="Power of 2 for grid size (size = 2^n + 1, default: 9 = 513x513)"
    )
    parser.add_argument(
        "-r", "--roughness", type=float, default=0.5,
        help="Roughness/decay factor (0.0-1.0, default: 0.5)"
    )
    parser.add_argument(
        "-o", "--output", type=str, default="terrain_ds.png",
        help="Output filename (default: terrain_ds.png)"
    )
    parser.add_argument("--seed", type=int, default=None, help="Random seed for reproducibility")
    parser.add_argument(
        "--grid", type=str, default=None,
        help="Generate stitched tile grid (e.g., '3x3'). True to arguments—no automatic adjustments."
    )
    parser.add_argument(
        "--grid2", type=str, default=None,
        help="Like --grid, but 'best effort' mode: adds corner smoothing, water border enforcement (with --border), and terrain stacking to achieve target ratios."
    )
    parser.add_argument(
        "--border", action="store_true",
        help="Pre-seed borders with low values (creates island-like terrain)"
    )
    parser.add_argument(
        "--water-ratio", type=float, default=None,
        help="Target water-to-land ratio (0.0-1.0, e.g., 0.3 for 30%% water)"
    )
    parser.add_argument(
        "--water-method", type=str, choices=["gamma", "piecewise"], default="gamma",
        help="Method for water ratio rescaling (default: gamma)"
    )
    parser.add_argument("--water-blur", action="store_true", help="Apply gaussian blur to water regions")
    parser.add_argument(
        "--smart-crop", action="store_true",
        help="Crop to minimize water border, then resize to original size"
    )
    
    args = parser.parse_args()
    size = (1 << args.power) + 1
    
    # Validate water ratio
    if args.water_ratio is not None and not (0.0 < args.water_ratio < 1.0):
        print(f"Error: --water-ratio must be between 0 and 1 (exclusive), got {args.water_ratio}")
        return
    
    grid_spec = args.grid or args.grid2
    if grid_spec:
        dims = parse_grid_spec(grid_spec)
        if dims is None:
            print(f"Error: Invalid grid format '{grid_spec}'. Use format like '3x3'")
            return
        grid_w, grid_h = dims
        
        mode_info = [f"{grid_w}x{grid_h} grid of {size}x{size} tiles", f"roughness={args.roughness}"]
        if args.border:
            mode_info.append("border seeding")
        if args.grid2:
            mode_info.append("grid2 mode")
        if args.water_ratio:
            mode_info.append(f"target water={args.water_ratio:.0%}")
        print(f"Generating {', '.join(mode_info)}")
        
        terrain = generate_grid_terrain(args, grid_w, grid_h)
        
        # Enforce water borders for --grid2 --border mode
        if args.grid2 and args.border:
            original_ratio = args.water_ratio or np.mean(terrain < WATER_THRESHOLD)
            terrain, success, final_ratio = enforce_water_border(terrain, args.water_method, original_ratio)
            if success and final_ratio != original_ratio:
                print(f"Border enforcement: Adjusted water ratio from {original_ratio:.3f} to {final_ratio:.3f} ({final_ratio - original_ratio:+.3f})")
            elif not success:
                print("Warning: Could not enforce water borders even at maximum water ratio")
            
            # Stack terrains if water ratio is too far from target
            if args.water_ratio is not None:
                current_ratio = np.mean(terrain < WATER_THRESHOLD)
                if current_ratio - args.water_ratio > STACKING_WATER_RATIO_TOLERANCE:
                    print(f"Water ratio {current_ratio:.3f} exceeds target {args.water_ratio:.3f} + tolerance {STACKING_WATER_RATIO_TOLERANCE:.0%}, stacking terrains...")
                    terrain, stack_count = stack_terrains_for_water_ratio(
                        terrain, args.water_ratio, args, grid_w, grid_h
                    )
                    if stack_count > 0:
                        final_ratio = np.mean(terrain < WATER_THRESHOLD)
                        print(f"Stacking complete: {stack_count} iteration(s), final water ratio: {final_ratio:.3f}")
        
        print(f"Total size: {terrain.shape[1]}x{terrain.shape[0]}")
    else:
        mode_info = [f"{size}x{size} terrain", f"roughness={args.roughness}"]
        if args.border:
            mode_info.append("border seeding")
        if args.water_ratio:
            mode_info.append(f"target water={args.water_ratio:.0%}")
        print(f"Generating {', '.join(mode_info)}")
        
        terrain = generate_single_terrain(args)
    
    # Report water ratio
    actual_ratio = np.mean(terrain < WATER_THRESHOLD)
    if args.water_ratio is not None:
        error = abs(actual_ratio - args.water_ratio)
        print(f"Water ratio - Target: {args.water_ratio:.3f}, Actual: {actual_ratio:.3f}, Error: {error:.4f}")
    else:
        print(f"Water ratio (natural): {actual_ratio:.3f}")
    
    # Apply postprocessing
    terrain = apply_postprocessing(terrain, args)
    
    # Re-enforce water borders after smart-crop if --grid2 --border was used
    if args.smart_crop and args.grid2 and args.border:
        original_ratio = np.mean(terrain < WATER_THRESHOLD)
        terrain, success, final_ratio = enforce_water_border(terrain, args.water_method, original_ratio)
        if success and final_ratio != original_ratio:
            print(f"Post-crop border enforcement: {original_ratio:.3f} -> {final_ratio:.3f} ({final_ratio - original_ratio:+.3f})")
        elif not success:
            print("Warning: Could not enforce water borders even at maximum water ratio")
    
    save_terrain_image(terrain, args.output)
    print("Done!")


 if __name__ == "__main__":
    main()
diff --git a/info.md b/info.md
Option	Description
`-n, --power`	Power of 2 for grid size (size = 2^n + 1, default: 9)
`-r, --roughness`	Roughness factor 0.0–1.0 (default: 0.5)
`-o, --output`	Output filename (default: terrain_ds.png)
`--seed`	Random seed for reproducibility
`--grid`	Generate stitched tile grid (e.g., "3x3"). True to arguments—no automatic adjustments.
`--grid2`	Like `--grid`, but "best effort" mode: adds corner smoothing, water border enforcement, and terrain stacking to achieve target ratios.
`--border`	Pre-seed borders with low values (island mode)
`--water-ratio`	Target water proportion 0.0–1.0
`--water-method`	Water rescaling: "gamma" (default) or "piecewise"
`--water-blur`	Apply gaussian blur to water regions (smooths coastlines)
`--smart-crop`	Crop to minimize water border, then resize to original size
No results found