export module ra3.terrain;

import std;

export import ra3.core;
export import ra3.fs;
export import ra3.render;
export import ra3.models;
import ra3.assets;

/**
 * The map's real terrain, read from the compiled `CkMp` chunk tree.
 *
 * RA3 maps carry a `HeightMapData` chunk (the elevation grid) and a
 * `BlendTileData` chunk (the per-cell tile index and the tile texture set).
 * This module parses both and rasterises the terrain from the user's own
 * terrain textures in `Data\Terrain.big` (RefPack + TGA), so a rendered map is
 * the actual in-game terrain rather than the `<map>_art.tga` overview.
 *
 * Format references (GPLv3): OpenSAGE `Data/Map/{HeightMapData,BlendTileData,
 * BlendTileTexture}.cs`, cross-checked against the retail maps.
 *
 * The RA3 tile value is `(cellIndex << 2) | variant`: the low two bits are the
 * blend variant, the rest indexes the map's global texture-cell table
 * (`TextureCellCount` cells, each texture owning `cellSize^2` of them).
 */
export namespace ra3::terrain {
    using ra3::core::uint8;
    using ra3::core::uint16;
    using ra3::core::uint32;
    using ra3::core::usize;
    using ra3::render::argb;
    using ra3::render::image;

    /** World units per terrain cell (SAGE `TheTerrainLogic`, `cellSize == 10`). */
    inline constexpr float cell_size = 10.0F;

    /** One entry of `BlendTileData`'s texture table. */
    struct texture_ref {
        uint32 cell_start = 0;
        uint32 cell_count = 0;
        uint32 cell_size = 0; ///< Cells per side; the source texture is `cell_size * 64` px.
        std::string name;
    };

    /**
     * One `BlendDescription` (`BlendTileData` tail): how a cell cross-fades into
     * a neighbour tile.
     *
     * `secondary_tile` is a packed tile value `(cellIndex << 2) | variant`, so
     * its cell is `secondary_tile >> 2`. `direction` is the SAGE `BlendDirection`
     * bitmask (1 right, 2 top, 4 top-right, 8 top-left); `flags` bit 0 flips the
     * direction and bit 1 marks a two-sided diagonal (see `shaders/terrain.frag`).
     */
    struct blend_description {
        uint16 secondary_tile = 0;
        uint8 direction = 0;
        uint8 flags = 0;
        bool two_sided = false;
    };

    /** The parsed terrain: elevation grid, tile grid, blends and tile textures. */
    struct map_data {
        uint32 width = 0;
        uint32 height = 0;
        uint32 border_width = 0;
        std::vector<uint16> elevations;
        std::vector<uint16> tiles;
        std::vector<uint16> blends;           ///< 1-based index into `blend_descriptions`, 0 = none.
        std::vector<uint16> three_way_blends; ///< Same, for the three-way corner blend.
        std::vector<blend_description> blend_descriptions;
        std::vector<texture_ref> textures;
        uint32 texture_cell_count = 0;
        bool has_water = false;
        float water_plane_z = 0.0F; ///< GlobalWaterSettings ReflectionPlaneZ, world units.
        bool valid = false;

        [[nodiscard]] auto elevation(uint32 x, uint32 y) const -> uint16 { return elevations[static_cast<usize>(y) * width + x]; }
        [[nodiscard]] auto tile(uint32 x, uint32 y) const -> uint16 { return tiles[static_cast<usize>(y) * width + x]; }
        [[nodiscard]] auto blend(uint32 x, uint32 y) const -> uint16 { return blends[static_cast<usize>(y) * width + x]; }
        [[nodiscard]] auto three_way_blend(uint32 x, uint32 y) const -> uint16 { return three_way_blends[static_cast<usize>(y) * width + x]; }

        /** World size in engine units. */
        [[nodiscard]] auto world_width() const -> float { return static_cast<float>(width) * cell_size; }
        [[nodiscard]] auto world_height() const -> float { return static_cast<float>(height) * cell_size; }
    };

    /** Thrown when the map's terrain chunks are malformed. */
    class terrain_error : public std::runtime_error {
    public:
        using std::runtime_error::runtime_error;
    };

    namespace detail {
        /** View a `uint8` range as bytes, the currency of `libra3assets`. */
        [[nodiscard]] inline auto as_bytes(std::span<const uint8> data) -> std::span<const std::byte> { return std::as_bytes(data); }

        /** View a `libra3assets` byte range as OpenRA3's `uint8`. */
        [[nodiscard]] inline auto as_u8(std::span<const std::byte> data) -> std::span<const uint8> {
            return {reinterpret_cast<const uint8 *>(data.data()), data.size()};
        }

        [[nodiscard]] inline auto read_u16(const uint8 *p) -> uint16 { return static_cast<uint16>(p[0]) | (static_cast<uint16>(p[1]) << 8U); }
        [[nodiscard]] inline auto read_u32(const uint8 *p) -> uint32 {
            return static_cast<uint32>(p[0]) | (static_cast<uint32>(p[1]) << 8U) | (static_cast<uint32>(p[2]) << 16U) |
                   (static_cast<uint32>(p[3]) << 24U);
        }

        // The `CkMp` container and `HeightMapData` chunk are modelled by
        // `libra3assets` (`ra3.assets`); `parse_map` reads them through a
        // `map_document`. Only `BlendTileData` (which the library does not type
        // yet) is parsed here, over the chunk payload the document exposes.

        /** Where the texture table ends and how many `BlendDescription`s follow. */
        struct texture_table_info {
            usize end = 0;           ///< Offset just past the last texture entry.
            uint32 blends_count = 0; ///< Raw `BlendDescriptions` count (0 or length+1).
        };

        /** Scan the tail of `BlendTileData` for the texture table. */
        [[nodiscard]] inline auto parse_textures(std::span<const uint8> data, usize from, usize end, map_data &out) -> texture_table_info {
            const auto printable = [](std::span<const uint8> s) {
                for (auto ch: s) {
                    if (ch < 32U || ch >= 127U) return false;
                }
                return true;
            };
            const auto valid = [&](usize q, texture_ref &ref, usize &next) -> bool {
                if (q + 18U > end) return false;
                const auto cell_start = read_u32(data.data() + q);
                const auto cell_count = read_u32(data.data() + q + 4U);
                const auto cell_side = read_u32(data.data() + q + 8U);
                const auto magic = read_u32(data.data() + q + 12U);
                if (cell_side == 0U || cell_side > 64U || magic != 0U || cell_count != cell_side * cell_side) return false;
                const auto name_len = read_u16(data.data() + q + 16U);
                if (name_len < 3U || name_len > 48U || q + 18U + name_len > end) return false;
                if (!printable(data.subspan(q + 18U, name_len))) return false;
                ref = {cell_start, cell_count, cell_side, std::string{reinterpret_cast<const char *>(data.data() + q + 18U), name_len}};
                next = q + 18U + name_len;
                return true;
            };

            for (usize q = from; q + 4U < end; ++q) {
                if (q < 12U) continue;
                const auto count = read_u32(data.data() + q);
                if (count == 0U || count > 256U) continue;
                usize cursor = q + 4U;
                std::vector<texture_ref> refs;
                bool ok = true;
                for (uint32 i = 0; i < count; ++i) {
                    texture_ref ref;
                    usize next = 0;
                    if (!valid(cursor, ref, next)) {
                        ok = false;
                        break;
                    }
                    refs.push_back(std::move(ref));
                    cursor = next;
                }
                if (!ok) continue;
                out.texture_cell_count = read_u32(data.data() + q - 12U);
                out.textures = std::move(refs);
                return texture_table_info{cursor, read_u32(data.data() + q - 8U)};
            }
            throw terrain_error("BlendTileData texture table not found");
        }

        /** Parse the `BlendDescriptions` that follow the texture table. */
        inline auto parse_blend_descriptions(std::span<const uint8> data, const texture_table_info &table, usize chunk_end, map_data &out) -> void {
            if (table.blends_count == 0U) return;
            const auto count = static_cast<uint32>(table.blends_count - 1U);
            usize p = table.end;
            if (p + 8U > chunk_end) return;
            p += 8U; // MagicValue1, MagicValue2
            out.blend_descriptions.reserve(count);
            for (uint32 i = 0; i < count && p + 18U <= chunk_end; ++i) {
                blend_description desc;
                desc.secondary_tile = static_cast<uint16>(read_u32(data.data() + p));
                const uint8 raw[4] = {data[p + 4U], data[p + 5U], data[p + 6U], data[p + 7U]};
                desc.flags = data[p + 8U];
                desc.two_sided = data[p + 9U] != 0U;
                uint8 direction = 0;
                for (uint32 k = 0; k < 4U; ++k) {
                    if (raw[k] != 0U) direction = static_cast<uint8>(direction | static_cast<uint8>(raw[k] << k));
                }
                desc.direction = direction;
                out.blend_descriptions.push_back(desc);
                p += 18U;
            }
        }

        inline auto parse_blend(std::span<const uint8> payload, uint16 version, map_data &out) -> void {
            const auto area = static_cast<usize>(out.width) * out.height;
            usize p = 0;
            const auto num_tiles = read_u32(payload.data() + p);
            p += 4;
            if (num_tiles != area) throw terrain_error("BlendTileData tile count mismatch");
            out.tiles.resize(area);
            for (usize i = 0; i < area; ++i) out.tiles[i] = read_u16(payload.data() + p + i * 2U);
            p += area * 2U;

            const auto bits = (version >= 14U && version < 24U) ? 32U : 16U;
            const auto word = bits / 8U;
            const auto read_index = [&](usize off) -> uint16 {
                return word == 4U ? static_cast<uint16>(read_u32(payload.data() + off)) : read_u16(payload.data() + off);
            };
            out.blends.resize(area);
            for (usize i = 0; i < area; ++i) out.blends[i] = read_index(p + i * word);
            p += area * word;
            out.three_way_blends.resize(area);
            for (usize i = 0; i < area; ++i) out.three_way_blends[i] = read_index(p + i * word);
            p += area * word;
            p += area * word; // CliffTextures (not rendered yet)

            const auto chunk_end = payload.size();
            const auto table = parse_textures(payload, p, chunk_end, out);
            parse_blend_descriptions(payload, table, chunk_end, out);
        }
    }

    /**
     * Parse the terrain chunks out of a `CkMp` (uncompressed) map payload.
     *
     * The container and `HeightMapData` come from `libra3assets`
     * (`map_document`); `BlendTileData` (which the library does not model yet)
     * is decoded here from the chunk payload the document exposes.
     *
     * @throws terrain_error if the chunk tree or terrain chunks are malformed.
     */
    [[nodiscard]] inline auto parse_map(std::span<const uint8> ckmp) -> map_data {
        try {
            auto document = ra3::assets::map_document::parse(detail::as_bytes(ckmp));

            map_data out;
            const auto height = document.height_map();
            if (!height) throw terrain_error("map has no HeightMapData chunk");
            out.width = height->width;
            out.height = height->height;
            out.border_width = height->border_width;
            out.elevations.assign(height->elevations.begin(), height->elevations.end());

            const auto *blend = document.find_chunk("BlendTileData");
            if (blend == nullptr) throw terrain_error("map has no BlendTileData chunk");
            detail::parse_blend(detail::as_u8(blend->payload), blend->version, out);

            if (const auto *water = document.find_chunk("GlobalWaterSettings"); water != nullptr && water->payload.size() >= 8U) {
                const auto *bytes = reinterpret_cast<const uint8 *>(water->payload.data());
                out.has_water = detail::read_u32(bytes) != 0U;
                const auto bits = detail::read_u32(bytes + 4U);
                std::memcpy(&out.water_plane_z, &bits, sizeof(out.water_plane_z));
            }
            out.valid = true;
            return out;
        } catch (const ra3::assets::asset_error &error) {
            throw terrain_error(error.what());
        }
    }

    /** The decoded terrain textures, parallel to `map_data::textures`. */
    struct texture_set {
        std::vector<image> images;

        [[nodiscard]] auto resolved() const -> usize {
            usize n = 0;
            for (const auto &img: images) {
                if (!img.empty()) ++n;
            }
            return n;
        }
    };

    namespace detail {
        /** Lowercased basename without `.tga`, e.g. `tdirt_yucatan03`. */
        [[nodiscard]] inline auto tga_stem(std::string_view path) -> std::string {
            const auto slash = path.find_last_of("\\/");
            auto name = std::string{path.substr(slash == std::string_view::npos ? 0U : slash + 1U)};
            if (name.size() > 4U) name.resize(name.size() - 4U);
            std::transform(name.begin(), name.end(), name.begin(), [](unsigned char ch) { return static_cast<char>(std::tolower(ch)); });
            return name;
        }
    }

    /**
     * Load the map's tile textures from the install's terrain archives.
     *
     * @param map The parsed map (its `textures` name the files).
     * @param data_dir The install's `Data` directory (holds `Terrain.big`).
     */
    [[nodiscard]] inline auto load_textures(const map_data &map, const std::filesystem::path &data_dir) -> texture_set {
        // Index every `art\terrain\*.tga` (minus normals) by lowercased stem.
        struct source {
            const fs::big_archive *archive = nullptr;
            std::string entry;
        };
        std::vector<fs::big_archive> archives;
        std::unordered_map<std::string, source> files;
        for (const auto &name: {"Terrain.big", "Core11.big"}) {
            const auto path = data_dir / name;
            std::error_code ec;
            if (!std::filesystem::exists(path, ec)) continue;
            archives.push_back(fs::big_archive::open(path));
        }
        for (const auto &archive: archives) {
            for (const auto &entry: archive.entries()) {
                auto stem = detail::tga_stem(entry.name);
                if (stem.size() > 4U && stem.ends_with("_nrm")) continue;
                if (!entry.name.ends_with(".tga") && !entry.name.ends_with(".TGA")) continue;
                files.try_emplace(stem, source{&archive, entry.name});
            }
        }

        texture_set set;
        set.images.resize(map.textures.size());
        for (usize i = 0; i < map.textures.size(); ++i) {
            auto name = map.textures[i].name;
            std::transform(name.begin(), name.end(), name.begin(), [](unsigned char ch) { return static_cast<char>(std::tolower(ch)); });
            const source *found = nullptr;
            for (const auto &candidate: {name, "t" + name, "tmisc_" + name}) {
                if (const auto it = files.find(candidate); it != files.end()) {
                    found = &it->second;
                    break;
                }
            }
            if (found == nullptr) {
                for (const auto &[stem, src]: files) {
                    if (stem.ends_with(name)) {
                        found = &src;
                        break;
                    }
                }
            }
            if (found == nullptr) continue;
            try {
                set.images[i] = ra3::render::decode_tga(found->archive->read(found->entry, true));
            } catch (const std::exception &) {
                // Leave the slot empty; the renderer falls back to a palette.
            }
        }
        return set;
    }

    /** Rasteriser parameters. */
    struct render_options {
        uint32 scale = 4;         ///< Output pixels per cell (horizontal).
        float pitch = 0.62F;      ///< Vertical foreshortening; 1 is straight down.
        float ambient = 0.40F;    ///< Ambient light fraction.
        std::array<float, 3> sun_dir{0.45F, 0.35F, 0.82F}; ///< From the scene toward the sun.
        float z_scale = 0.0390625F; ///< HeightMapData vertical scale (version >= 5).
        uint32 water_color = argb(34, 74, 118);
        float water_alpha = 0.60F;
        float shadow_strength = 0.45F; ///< Fraction of light removed where a shadow falls.
        uint32 shadow_steps = 24U;
        uint32 supersample = 2U;       ///< Render at Nx and box-downsample (antialiasing).
    };

    namespace detail {
        /** Index `*.tga` under a terrain dir by stem (lower-cased), ignoring normals. */
        [[nodiscard]] inline auto terrain_file_index(const std::filesystem::path &dir)
            -> std::unordered_map<std::string, std::filesystem::path> {
            std::unordered_map<std::string, std::filesystem::path> files;
            std::error_code ec;
            if (!std::filesystem::is_directory(dir, ec)) return files;
            // Recursive: our own extract writes flat `terrain/*.tga`, ra3tools
            // writes nested `.../art/terrain/*.tga`. Only index `terrain` dirs
            // so a full asset dump does not pull in every unrelated TGA.
            for (const auto &entry: std::filesystem::recursive_directory_iterator(dir, ec)) {
                if (!entry.is_regular_file()) continue;
                const auto &path = entry.path();
                auto parent = path.parent_path().filename().string();
                std::transform(parent.begin(), parent.end(), parent.begin(), [](unsigned char ch) { return static_cast<char>(std::tolower(ch)); });
                if (parent != "terrain") continue;
                auto stem = tga_stem(path.filename().string());
                if (stem.size() > 4U && stem.ends_with("_nrm")) continue;
                files.try_emplace(stem, path);
            }
            return files;
        }

        /** Match one map texture name against the index (see `load_textures_from_dir`). */
        [[nodiscard]] inline auto match_terrain_file(const std::unordered_map<std::string, std::filesystem::path> &files, std::string name)
            -> std::filesystem::path {
            std::transform(name.begin(), name.end(), name.begin(), [](unsigned char ch) { return static_cast<char>(std::tolower(ch)); });
            for (const auto &candidate: {name, "t" + name, "tmisc_" + name}) {
                if (const auto it = files.find(candidate); it != files.end()) return it->second;
            }
            for (const auto &[stem, path]: files) {
                if (stem.ends_with(name)) return path;
            }
            return {};
        }
    } // namespace detail

    /**
     * The loose `*.tga` files `map` resolves to under `dir`, de-duplicated.
     *
     * Used to stage just the tiles a single map needs (e.g. the wasm preload).
     */
    [[nodiscard]] inline auto resolve_texture_files(const map_data &map, const std::filesystem::path &dir) -> std::vector<std::filesystem::path> {
        const auto files = detail::terrain_file_index(dir);
        std::vector<std::filesystem::path> resolved;
        for (const auto &texture: map.textures) {
            if (auto found = detail::match_terrain_file(files, texture.name); !found.empty()) resolved.push_back(std::move(found));
        }
        std::sort(resolved.begin(), resolved.end());
        resolved.erase(std::unique(resolved.begin(), resolved.end()), resolved.end());
        return resolved;
    }

    /** Load tile textures from a directory of loose `*.tga` files (extracted assets). */
    [[nodiscard]] inline auto load_textures_from_dir(const map_data &map, const std::filesystem::path &dir,
                                                    const std::function<void(float)> &progress = {}) -> texture_set {
        const auto files = detail::terrain_file_index(dir);
        texture_set set;
        set.images.resize(map.textures.size());
        for (usize i = 0; i < map.textures.size(); ++i) {
            const auto found = detail::match_terrain_file(files, map.textures[i].name);
            if (!found.empty()) {
                try {
                    std::ifstream in(found, std::ios::binary);
                    std::vector<uint8> raw((std::istreambuf_iterator<char>(in)), std::istreambuf_iterator<char>());
                    set.images[i] = ra3::render::decode_tga(raw);
                } catch (const std::exception &) {
                }
            }
            if (progress && !map.textures.empty()) progress(static_cast<float>(i + 1U) / static_cast<float>(map.textures.size()));
        }
        return set;
    }

    /** Terrain packed for GPU upload: heights, per-cell blend data, texture array. */
    struct gpu_terrain {
        uint32 width = 0;
        uint32 height = 0;
        std::vector<uint16> heights;   ///< W*H raw elevations (R16).
        std::vector<uint16> cell_data; ///< W*H * 4 u16: base layer, blend layer, three-way layer, packed blend info.
        std::vector<uint32> layers;    ///< `layer_count` square RGBA8 textures laid out back to back.
        uint32 layer_count = 0;
        uint32 layer_size = 0;
        float cell_span = 8.0F; ///< Cells per texture repeat (`2 * cellSize`), for `uv = cell / span`.
        float z_scale = 0.0390625F;
        bool has_water = false;
        float water_z = 0.0F;
        ra3::models::scene objects; ///< Buildings and props placed on the map (world-space triangle soup).
    };

    /**
     * Build the GPU textures for `map` (heights, per-cell blend data, texture array).
     *
     * The tile textures are uploaded into an array and sampled **continuously**
     * (`uv = cell / (2 * cellSize)`, SAGE/OpenSAGE `Terrain.frag`), so the
     * texture never restarts at a cell edge. The per-cell record carries the
     * base texture layer, the blend (and three-way) secondary layer and the
     * packed blend direction/flags, which the shader ramps across the cell.
     */
    [[nodiscard]] inline auto build_gpu_terrain(const map_data &map, const texture_set &set, const render_options &options,
                                                const ra3::models::scene &objects, const std::function<void(float)> &progress = {}) -> gpu_terrain {
        gpu_terrain out;
        out.width = map.width;
        out.height = map.height;
        out.z_scale = options.z_scale;
        out.has_water = map.has_water;
        out.water_z = map.water_plane_z;
        out.objects = objects;

        // Cell index -> texture layer.
        const auto layer_of = [&](uint32 cell_index) -> uint16 {
            for (usize i = 0; i < map.textures.size(); ++i) {
                const auto &t = map.textures[i];
                if (cell_index >= t.cell_start && cell_index < t.cell_start + t.cell_count) return static_cast<uint16>(i);
            }
            return 0xFFFFU;
        };

        const auto area = static_cast<usize>(map.width) * map.height;
        out.heights.resize(area);
        out.cell_data.assign(area * 4U, 0U);
        for (usize i = 0; i < area; ++i) {
            out.heights[i] = map.elevations[i];
            const auto base = layer_of(map.tiles[i] >> 2U);
            auto blend_layer = base;
            auto three_layer = base;
            uint32 packed = 0;
            const auto resolve = [&](uint16 index, uint32 shift) -> uint16 {
                if (index == 0U || index > map.blend_descriptions.size()) return base;
                const auto &desc = map.blend_descriptions[index - 1U];
                const auto layer = layer_of(desc.secondary_tile >> 2U);
                uint32 flags = desc.flags & 0x1U;
                if (desc.two_sided) flags |= 0x2U;
                packed |= static_cast<uint32>(static_cast<uint16>((desc.direction & 0xFU) | ((flags & 0x3U) << 4U))) << shift;
                return layer;
            };
            if (i < map.blends.size()) blend_layer = resolve(map.blends[i], 0U);
            if (i < map.three_way_blends.size()) three_layer = resolve(map.three_way_blends[i], 8U);
            out.cell_data[i * 4U + 0U] = base;
            out.cell_data[i * 4U + 1U] = blend_layer;
            out.cell_data[i * 4U + 2U] = three_layer;
            out.cell_data[i * 4U + 3U] = static_cast<uint16>(packed);
        }

        out.layer_count = static_cast<uint32>(std::max<usize>(1U, map.textures.size()));
        uint32 layer_size = 64U;
        for (const auto &img: set.images) {
            if (!img.empty()) layer_size = std::max(layer_size, img.width());
        }
        out.layer_size = layer_size;
        out.cell_span = 8.0F;
        for (const auto &ref: map.textures) {
            if (ref.cell_size != 0U) {
                out.cell_span = 2.0F * static_cast<float>(ref.cell_size);
                break;
            }
        }

        out.layers.assign(static_cast<usize>(out.layer_count) * layer_size * layer_size, 0xFF3A4550U);
        for (usize i = 0; i < map.textures.size(); ++i) {
            const auto &img = set.images[i];
            if (img.empty()) continue;
            for (uint32 y = 0; y < layer_size; ++y) {
                const auto sy = std::min(img.height() - 1U, y * img.height() / layer_size);
                for (uint32 x = 0; x < layer_size; ++x) {
                    const auto sx = std::min(img.width() - 1U, x * img.width() / layer_size);
                    out.layers[(i * layer_size + y) * layer_size + x] = img.data()[static_cast<usize>(sy) * img.width() + sx];
                }
            }
        }
        if (progress) progress(1.0F);
        return out;
    }

    /** Terrain without any placed objects. */
    [[nodiscard]] inline auto build_gpu_terrain(const map_data &map, const texture_set &set, const render_options &options = {},
                                               const std::function<void(float)> &progress = {}) -> gpu_terrain {
        return build_gpu_terrain(map, set, options, ra3::models::scene{}, progress);
    }

    namespace detail {
        /**
         * The source texture a tile cell maps to. The texture is sampled
         * **continuously** (SAGE `uv / (cellSize * 2)`), so adjacent cells of
         * the same material do not restart the texture — that restart is what
         * left a seam along every cell edge.
         */
        struct layer_ref {
            const image *img = nullptr;
            float span = 8.0F; ///< Cells per texture repeat (`2 * cellSize`).
            uint32 fallback = 0;
        };

        [[nodiscard]] inline auto resolve_layer(const map_data &map, const texture_set &set, uint32 cell_index) -> layer_ref {
            for (usize i = 0; i < map.textures.size(); ++i) {
                const auto &t = map.textures[i];
                if (cell_index < t.cell_start || cell_index >= t.cell_start + t.cell_count) continue;
                const auto &img = set.images[i];
                const auto seed = static_cast<uint32>(i);
                if (img.empty() || t.cell_size == 0U) return {nullptr, 8.0F, seed};
                return {&img, 2.0F * static_cast<float>(t.cell_size), seed};
            }
            return {nullptr, 8.0F, cell_index};
        }

        [[nodiscard]] inline auto palette(uint32 seed) -> uint32 {
            const auto h = seed * 2654435761U;
            return argb(static_cast<uint8>(70U + (h & 0x5FU)), static_cast<uint8>(70U + ((h >> 8U) & 0x5FU)),
                        static_cast<uint8>(70U + ((h >> 16U) & 0x5FU)));
        }

        /**
         * SAGE blend ramp. Returns 0 on one edge of the cell and 1 on the
         * opposite edge, so the base tile cross-fades into the secondary tile
         * across the whole cell rather than at a hard line. Mirrors
         * `CalculateBlendFactor` in the retail `Terrain.fx` / OpenSAGE
         * `Terrain.frag`.
         */
        [[nodiscard]] inline auto blend_factor(uint32 direction, uint32 flags, float fx, float fy) -> float {
            const bool flipped = (flags & 1U) != 0U;
            const bool two_sided = (flags & 2U) != 0U;
            if (flipped) {
                if (direction == 1U) {
                    fx = 1.0F - fx;
                } else if (direction == 2U || direction == 4U || direction == 8U) {
                    fy = 1.0F - fy;
                }
            }
            const auto sat = [](float v) { return std::clamp(v, 0.0F, 1.0F); };
            switch (direction) {
            case 1U: return fx;
            case 2U: return fy;
            case 4U: {
                const auto s = (1.0F - fx) + (1.0F - fy);
                return two_sided ? 1.0F - sat(s - 1.0F) : sat(1.0F - s);
            }
            case 8U: {
                const auto s = fx + (1.0F - fy);
                return two_sided ? 1.0F - sat(s - 1.0F) : sat(1.0F - s);
            }
            default: return 0.0F;
            }
        }

        /** Pack `(direction, flags)` for the shader (4-bit direction + 2 flag bits). */
        [[nodiscard]] inline auto pack_blend(uint32 direction, uint32 flags) -> uint16 {
            return static_cast<uint16>((direction & 0xFU) | ((flags & 0x3U) << 4U));
        }

        [[nodiscard]] inline auto wrap01(float value) -> float { return value - std::floor(value); }

        /** Nearest sample of a layer at global cell coordinates `(wx, wy)`. */
        [[nodiscard]] inline auto sample_layer(const layer_ref &layer, float wx, float wy) -> uint32 {
            if (layer.img == nullptr || layer.img->empty()) return palette(layer.fallback);
            const auto iw = layer.img->width();
            const auto ih = layer.img->height();
            const auto sx = std::min(iw - 1U, static_cast<uint32>(wrap01(wx / layer.span) * static_cast<float>(iw)));
            const auto sy = std::min(ih - 1U, static_cast<uint32>(wrap01(wy / layer.span) * static_cast<float>(ih)));
            return layer.img->data()[static_cast<usize>(sy) * iw + sx];
        }

        /** Bilinear sample of a layer at global cell coordinates `(wx, wy)`. */
        [[nodiscard]] inline auto sample_layer_bilinear(const layer_ref &layer, float wx, float wy) -> uint32 {
            if (layer.img == nullptr || layer.img->empty()) return palette(layer.fallback);
            const auto iw = layer.img->width();
            const auto ih = layer.img->height();
            const auto px = static_cast<float>(iw);
            const auto py = static_cast<float>(ih);
            const auto sx = wrap01(wx / layer.span) * px - 0.5F;
            const auto sy = wrap01(wy / layer.span) * py - 0.5F;
            const auto x0 = static_cast<int>(std::floor(sx));
            const auto y0 = static_cast<int>(std::floor(sy));
            const auto tx = sx - static_cast<float>(x0);
            const auto ty = sy - static_cast<float>(y0);
            const auto wrap_idx = [](int v, uint32 n) -> uint32 { return static_cast<uint32>(((v % static_cast<int>(n)) + static_cast<int>(n)) % static_cast<int>(n)); };
            const auto texel = [&](int X, int Y) -> uint32 {
                return layer.img->data()[static_cast<usize>(wrap_idx(Y, ih)) * iw + wrap_idx(X, iw)];
            };
            const auto c00 = texel(x0, y0);
            const auto c10 = texel(x0 + 1, y0);
            const auto c01 = texel(x0, y0 + 1);
            const auto c11 = texel(x0 + 1, y0 + 1);
            const auto mix = [](float a, float b, float t) { return a + (b - a) * t; };
            const auto top = [&](uint32 sft) { return mix(static_cast<float>((c00 >> sft) & 0xFFU), static_cast<float>((c10 >> sft) & 0xFFU), tx); };
            const auto bot = [&](uint32 sft) { return mix(static_cast<float>((c01 >> sft) & 0xFFU), static_cast<float>((c11 >> sft) & 0xFFU), tx); };
            return argb(static_cast<uint8>(mix(top(16), bot(16), ty)), static_cast<uint8>(mix(top(8), bot(8), ty)),
                        static_cast<uint8>(mix(top(0), bot(0), ty)));
        }

        [[nodiscard]] inline auto mix_color(uint32 a, uint32 b, float w) -> uint32 {
            const auto mix = [w](uint32 x, uint32 y, uint32 shift) {
                return static_cast<uint32>(static_cast<float>((x >> shift) & 0xFFU) * (1.0F - w) + static_cast<float>((y >> shift) & 0xFFU) * w);
            };
            return argb(static_cast<uint8>(mix(a, b, 16)), static_cast<uint8>(mix(a, b, 8)), static_cast<uint8>(mix(a, b, 0)));
        }

        /** The three layers a map cell can draw: base, blend and three-way blend. */
        struct cell_layers {
            layer_ref base;
            layer_ref blend;
            layer_ref three;
            uint16 packed = 0; ///< Same layout as `gpu_terrain::cell_data[.w]`.
        };

        [[nodiscard]] inline auto resolve_layers(const map_data &map, const texture_set &set, usize index) -> cell_layers {
            cell_layers out;
            out.base = resolve_layer(map, set, static_cast<uint32>(map.tiles[index] >> 2U));
            out.blend = out.base;
            out.three = out.base;
            const auto apply = [&](uint16 blend_index, uint32 shift, layer_ref &target) {
                if (blend_index == 0U || blend_index > map.blend_descriptions.size()) return;
                const auto &desc = map.blend_descriptions[blend_index - 1U];
                target = resolve_layer(map, set, static_cast<uint32>(desc.secondary_tile >> 2U));
                uint32 flags = desc.flags & 0x1U;
                if (desc.two_sided) flags |= 0x2U;
                out.packed = static_cast<uint16>(out.packed | (pack_blend(desc.direction, flags) << shift));
            };
            if (index < map.blends.size()) apply(map.blends[index], 0U, out.blend);
            if (index < map.three_way_blends.size()) apply(map.three_way_blends[index], 8U, out.three);
            return out;
        }

        /** Blend a cell's layers at global `(wx, wy)` with the ramp at local `(fx, fy)`. */
        [[nodiscard]] inline auto blend_layers(const cell_layers &layers, float wx, float wy, float fx, float fy) -> uint32 {
            const auto dir1 = static_cast<uint32>(layers.packed & 0xFU);
            const auto flags1 = static_cast<uint32>((layers.packed >> 4U) & 0x3U);
            const auto dir2 = static_cast<uint32>((layers.packed >> 8U) & 0xFU);
            const auto flags2 = static_cast<uint32>((layers.packed >> 12U) & 0x3U);
            const auto f1 = blend_factor(dir1, flags1, fx, fy);
            const auto f2 = blend_factor(dir2, flags2, fx, fy);
            const auto c0 = sample_layer(layers.base, wx, wy);
            const auto c1 = f1 > 0.0F ? sample_layer(layers.blend, wx, wy) : c0;
            const auto c2 = f2 > 0.0F ? sample_layer(layers.three, wx, wy) : c0;
            return mix_color(mix_color(c0, c1, f1), c2, f2);
        }

        /**
         * Rasterise the map's building/prop scene over an already ray-marched
         * terrain image, depth-testing against it.
         *
         * The camera basis is the one `render3d` used, so the two passes agree;
         * `zbuf` holds the terrain's view-space depth per pixel (large where the
         * ray hit nothing). Triangles are z-tested and perspective-correct.
         */
        inline auto rasterize_objects(image &hi, std::vector<float> &zbuf, const ra3::models::scene &scene, const std::array<float, 3> &cam,
                                      const std::array<float, 3> &f, const std::array<float, 3> &r, const std::array<float, 3> &u, float tan_half,
                                      float aspect) -> void {
            const auto rw = static_cast<int>(hi.width());
            const auto rh = static_cast<int>(hi.height());
            if (rw <= 0 || rh <= 0) return;
            const auto dot3 = [](const std::array<float, 3> &a, const std::array<float, 3> &b) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; };
            float sun[3] = {0.45F, 0.35F, 0.82F};
            const auto sl = std::sqrt(sun[0] * sun[0] + sun[1] * sun[1] + sun[2] * sun[2]);
            sun[0] /= sl;
            sun[1] /= sl;
            sun[2] /= sl;
            constexpr float ambient = 0.38F;
            constexpr float near_plane = 10.0F;

            const auto sample = [&](uint32 layer, float tu, float tv) -> uint32 {
                if (layer >= scene.textures.size() || scene.textures[layer].empty()) return argb(140, 140, 140);
                const auto &img = scene.textures[layer];
                const auto wrap = [](float x) { return x - std::floor(x); };
                const auto sx = std::min(img.width() - 1U, static_cast<uint32>(wrap(tu) * static_cast<float>(img.width())));
                const auto sy = std::min(img.height() - 1U, static_cast<uint32>(wrap(tv) * static_cast<float>(img.height())));
                return img.data()[static_cast<usize>(sy) * img.width() + sx];
            };

            // Average colour per texture, used when a triangle covers fewer pixels
            // than texels (minification) — a cheap mip-0-away fallback that keeps
            // distant props from shimmering.
            std::vector<uint32> average(scene.textures.size(), argb(140, 140, 140));
            for (usize i = 0; i < scene.textures.size(); ++i) {
                const auto &img = scene.textures[i];
                if (img.empty()) continue;
                usize r = 0;
                usize g = 0;
                usize b = 0;
                for (usize p = 0; p < static_cast<usize>(img.width()) * img.height(); ++p) {
                    r += (img.data()[p] >> 16U) & 0xFFU;
                    g += (img.data()[p] >> 8U) & 0xFFU;
                    b += img.data()[p] & 0xFFU;
                }
                const auto n = static_cast<usize>(img.width()) * img.height();
                average[i] = argb(static_cast<uint8>(r / n), static_cast<uint8>(g / n), static_cast<uint8>(b / n));
            }

            struct projected {
                float sx = 0.0F;
                float sy = 0.0F;
                float inv_a = 0.0F; ///< 1 / view-space depth
            };
            const auto project = [&](const ra3::models::vertex &v, projected &out) -> bool {
                const std::array<float, 3> rel{v.x - cam[0], v.y - cam[1], v.z - cam[2]};
                const auto a = dot3(rel, f);
                if (a <= near_plane) return false;
                const auto ndc_x = (dot3(rel, r) / a) / (tan_half * aspect);
                const auto ndc_y = (dot3(rel, u) / a) / tan_half;
                out.sx = (ndc_x * 0.5F + 0.5F) * static_cast<float>(rw);
                out.sy = (0.5F - ndc_y * 0.5F) * static_cast<float>(rh);
                out.inv_a = 1.0F / a;
                return true;
            };

            for (usize t = 0; t + 2U < scene.indices.size(); t += 3U) {
                const auto &v0 = scene.vertices[scene.indices[t]];
                const auto &v1 = scene.vertices[scene.indices[t + 1U]];
                const auto &v2 = scene.vertices[scene.indices[t + 2U]];
                projected p0, p1, p2;
                if (!project(v0, p0) || !project(v1, p1) || !project(v2, p2)) continue;

                const auto area = (p1.sx - p0.sx) * (p2.sy - p0.sy) - (p1.sy - p0.sy) * (p2.sx - p0.sx);
                if (std::abs(area) < 1.0e-6F) continue;
                const auto sign = area < 0.0F ? -1.0F : 1.0F;

                // Texture footprint: if the triangle covers more texels than
                // pixels it is minified, so fall back to the texture average.
                const auto uv_area = std::abs((v1.u - v0.u) * (v2.v - v0.v) - (v2.u - v0.u) * (v1.v - v0.v));
                const auto texel_footprint = uv_area * static_cast<float>(scene.texture_size) * static_cast<float>(scene.texture_size);
                const auto minified = texel_footprint > 2.0F * std::abs(area);
                const auto flat_layer = static_cast<uint32>(v0.layer + 0.5F);
                const auto flat_color = flat_layer < average.size() ? average[flat_layer] : argb(140, 140, 140);

                const auto min_x = std::max(0, static_cast<int>(std::floor(std::min({p0.sx, p1.sx, p2.sx}))));
                const auto max_x = std::min(rw - 1, static_cast<int>(std::ceil(std::max({p0.sx, p1.sx, p2.sx}))));
                const auto min_y = std::max(0, static_cast<int>(std::floor(std::min({p0.sy, p1.sy, p2.sy}))));
                const auto max_y = std::min(rh - 1, static_cast<int>(std::ceil(std::max({p0.sy, p1.sy, p2.sy}))));

                for (int y = min_y; y <= max_y; ++y) {
                    for (int x = min_x; x <= max_x; ++x) {
                        const auto px = static_cast<float>(x) + 0.5F;
                        const auto py = static_cast<float>(y) + 0.5F;
                        auto w0 = ((p1.sx - p0.sx) * (py - p0.sy) - (p1.sy - p0.sy) * (px - p0.sx)) * sign;
                        auto w1 = ((p2.sx - p1.sx) * (py - p1.sy) - (p2.sy - p1.sy) * (px - p1.sx)) * sign;
                        auto w2 = ((p0.sx - p2.sx) * (py - p2.sy) - (p0.sy - p2.sy) * (px - p2.sx)) * sign;
                        if (w0 < 0.0F || w1 < 0.0F || w2 < 0.0F) continue;
                        const auto sum = w0 + w1 + w2;
                        if (sum <= 0.0F) continue;
                        w0 /= sum;
                        w1 /= sum;
                        w2 /= sum;

                        // Perspective-correct depth and attributes.
                        const auto inv_a = w0 * p0.inv_a + w1 * p1.inv_a + w2 * p2.inv_a;
                        const auto depth = 1.0F / inv_a;
                        const auto pixel = static_cast<usize>(y) * static_cast<usize>(rw) + static_cast<usize>(x);
                        // Ground decals carry a per-vertex bias toward the camera
                        // (the retail shader's screen-space depth offset); it
                        // keeps roads/sidewalks from z-fighting the terrain.
                        const auto bias = w0 * v0.bias + w1 * v1.bias + w2 * v2.bias;
                        if (depth - bias - 1.0F >= zbuf[pixel]) continue;

                        const auto tu = (w0 * v0.u * p0.inv_a + w1 * v1.u * p1.inv_a + w2 * v2.u * p2.inv_a) * depth;
                        const auto tv = (w0 * v0.v * p0.inv_a + w1 * v1.v * p1.inv_a + w2 * v2.v * p2.inv_a) * depth;
                        float nx = w0 * v0.nx + w1 * v1.nx + w2 * v2.nx;
                        float ny = w0 * v0.ny + w1 * v1.ny + w2 * v2.ny;
                        float nz = w0 * v0.nz + w1 * v1.nz + w2 * v2.nz;
                        const auto nl = std::sqrt(nx * nx + ny * ny + nz * nz);
                        if (nl > 1.0e-6F) {
                            nx /= nl;
                            ny /= nl;
                            nz /= nl;
                        }
                        const auto lambert = std::max(0.0F, std::abs(nx * sun[0] + ny * sun[1] + nz * sun[2]));
                        const auto shade = ambient + (1.0F - ambient) * lambert;
                        const auto texel = minified ? flat_color : sample(flat_layer, tu, tv);
                        const auto mod = [&](uint32 shift) {
                            return static_cast<uint8>(std::clamp(static_cast<float>((texel >> shift) & 0xFFU) * shade, 0.0F, 255.0F));
                        };
                        hi.data()[pixel] = argb(mod(16U), mod(8U), mod(0U));
                        zbuf[pixel] = depth;
                    }
                }
            }
        }
    }

    /**
     * Rasterise the terrain into an oblique, textured, lit image.
     *
     * Rows are compressed by `options.pitch` to emulate RA3's pitched tactical
     * view. Each cell is sampled from its 64 px tile texture (so the real map
     * detail shows rather than an averaged block); a directional sun shades the
     * relief (Lambert plus a heightfield shadow march); and cells below the
     * map's water plane (`GlobalWaterSettings` `ReflectionPlaneZ`) are covered
     * with translucent water.
     */
    [[nodiscard]] inline auto render(const map_data &map, const texture_set &set, const render_options &options = {}) -> image {
        if (!map.valid) throw terrain_error("terrain not parsed");
        const auto scale = std::max(1U, options.scale);
        const auto pitch = std::clamp(options.pitch, 0.05F, 1.0F);
        const auto out_w = map.width * scale;
        const auto out_h = static_cast<uint32>(std::lround(static_cast<double>(map.height) * scale * pitch));

        // Directional light + a cheap heightfield shadow.
        float lx = options.sun_dir[0];
        float ly = options.sun_dir[1];
        float lz = options.sun_dir[2];
        const auto len = std::sqrt(lx * lx + ly * ly + lz * lz);
        lx /= len;
        ly /= len;
        lz /= len;
        const auto horiz = std::max(1.0e-4F, std::sqrt(lx * lx + ly * ly));
        const auto sun_step_x = lx / horiz;
        const auto sun_step_y = ly / horiz;

        const auto z_of = [&](int x, int y) -> float {
            const auto cx = std::clamp(x, 0, static_cast<int>(map.width) - 1);
            const auto cy = std::clamp(y, 0, static_cast<int>(map.height) - 1);
            return static_cast<float>(map.elevation(static_cast<uint32>(cx), static_cast<uint32>(cy))) * options.z_scale;
        };

        std::vector<detail::cell_layers> layers(static_cast<usize>(map.width) * map.height);
        std::vector<float> light(static_cast<usize>(map.width) * map.height);
        for (uint32 y = 0; y < map.height; ++y) {
            for (uint32 x = 0; x < map.width; ++x) {
                const auto index = static_cast<usize>(y) * map.width + x;
                layers[index] = detail::resolve_layers(map, set, index);

                const auto nx = z_of(static_cast<int>(x) - 1, static_cast<int>(y)) - z_of(static_cast<int>(x) + 1, static_cast<int>(y));
                const auto ny = z_of(static_cast<int>(x), static_cast<int>(y) - 1) - z_of(static_cast<int>(x), static_cast<int>(y) + 1);
                const auto nlen = std::sqrt(nx * nx + ny * ny + 4.0F * cell_size * cell_size);
                const auto nz_norm = 2.0F * cell_size / nlen;
                auto lambert = (nx / nlen) * lx + (ny / nlen) * ly + nz_norm * lz;
                lambert = std::max(0.0F, lambert);

                // Shadow: march toward the sun; if terrain rises above the ray, shade.
                const auto z0 = z_of(static_cast<int>(x), static_cast<int>(y));
                for (uint32 s = 1; s <= options.shadow_steps; ++s) {
                    const auto sx = static_cast<int>(x) + static_cast<int>(std::lround(sun_step_x * static_cast<float>(s)));
                    const auto sy = static_cast<int>(y) + static_cast<int>(std::lround(sun_step_y * static_cast<float>(s)));
                    if (sx < 0 || sy < 0 || sx >= static_cast<int>(map.width) || sy >= static_cast<int>(map.height)) break;
                    const auto ray_z = z0 + (lz / horiz) * static_cast<float>(s) * cell_size;
                    if (z_of(sx, sy) > ray_z) {
                        lambert *= options.shadow_strength;
                        break;
                    }
                }

                auto shade = options.ambient + (1.0F - options.ambient) * lambert;
                const auto z_world = z0;
                if (map.has_water && z_world < map.water_plane_z) shade *= 0.75F; // dim the seabed
                light[index] = std::clamp(shade, 0.0F, 1.5F);
            }
        }

        image out(out_w, std::max(1U, out_h), argb(0, 0, 0));
        const auto *pixels = out.data();
        (void)pixels;

        for (uint32 oy = 0; oy < out_h; ++oy) {
            const auto fcy = static_cast<float>(oy) / (static_cast<float>(scale) * pitch);
            const auto cy = std::min(map.height - 1U, static_cast<uint32>(fcy));
            const auto fy = fcy - static_cast<float>(cy);
            for (uint32 ox = 0; ox < out_w; ++ox) {
                const auto fcx = static_cast<float>(ox) / static_cast<float>(scale);
                const auto cx = std::min(map.width - 1U, static_cast<uint32>(fcx));
                const auto fx = fcx - static_cast<float>(cx);
                const auto index = static_cast<usize>(cy) * map.width + cx;

                const auto base = detail::blend_layers(layers[index], fcx, fcy, fx, fy);

                const auto shade = light[index];
                const auto modulate = [&](uint32 c) { return static_cast<uint32>(std::clamp(static_cast<float>(c) * shade, 0.0F, 255.0F)); };
                uint32 r = modulate((base >> 16U) & 0xFFU);
                uint32 g = modulate((base >> 8U) & 0xFFU);
                uint32 b = modulate(base & 0xFFU);

                if (map.has_water && z_of(static_cast<int>(cx), static_cast<int>(cy)) < map.water_plane_z) {
                    const auto wr = (options.water_color >> 16U) & 0xFFU;
                    const auto wg = (options.water_color >> 8U) & 0xFFU;
                    const auto wb = options.water_color & 0xFFU;
                    const auto a = options.water_alpha;
                    r = static_cast<uint32>(r * (1.0F - a) + static_cast<float>(wr) * a);
                    g = static_cast<uint32>(g * (1.0F - a) + static_cast<float>(wg) * a);
                    b = static_cast<uint32>(b * (1.0F - a) + static_cast<float>(wb) * a);
                }

                const auto color = argb(static_cast<uint8>(r), static_cast<uint8>(g), static_cast<uint8>(b));
                out.data()[static_cast<usize>(oy) * out_w + ox] = color;
            }
        }
        return out;
    }

    /** Map a world position to a pixel in the terrain image (row 0 == max Y). */
    [[nodiscard]] inline auto world_to_pixel(const map_data &map, double world_x, double world_y, uint32 scale, float pitch) -> std::pair<int, int> {
        const auto cx = world_x / cell_size;
        const auto cy = static_cast<double>(map.height) - world_y / cell_size;
        return {static_cast<int>(std::lround(cx * scale)), static_cast<int>(std::lround(cy * scale * pitch))};
    }

    /** A perspective camera aimed at a ground target, like the retail tactical view. */
    using ra3::render::camera3d;

    /**
     * Render a true 3D view of the heightfield with a perspective camera.
     *
     * Rays are marched against the elevation grid; on a hit the tile's 64 px
     * texture is sampled at the exact hit fractional coordinate (so the map
     * stays sharp at any zoom), shaded by a directional sun with a shadow march,
     * and water cells get the water plane colour. This is real 3D geometry, not
     * a transformed 2D image.
     */
    [[nodiscard]] inline auto render3d(const map_data &map, const texture_set &set, const camera3d &camera, uint32 out_w, uint32 out_h,
                                       const render_options &options = {}, const ra3::models::scene *objects = nullptr) -> image {
        if (!map.valid) throw terrain_error("terrain not parsed");
        out_w = std::max(1U, out_w);
        out_h = std::max(1U, out_h);
        const auto ss = std::clamp(options.supersample, 1U, 4U);
        const auto rw = out_w * ss;
        const auto rh = out_h * ss;

        const auto world_w = static_cast<float>(map.width) * cell_size;
        const auto world_h = static_cast<float>(map.height) * cell_size;
        const auto zs = options.z_scale;

        std::vector<float> height(static_cast<usize>(map.width) * map.height);
        std::vector<detail::cell_layers> layers(height.size());
        for (uint32 y = 0; y < map.height; ++y) {
            for (uint32 x = 0; x < map.width; ++x) {
                const auto i = static_cast<usize>(y) * map.width + x;
                height[i] = static_cast<float>(map.elevation(x, y)) * zs;
                layers[i] = detail::resolve_layers(map, set, i);
            }
        }
        const auto cell_at = [&](int cx, int cy) -> const detail::cell_layers & {
            cx = std::clamp(cx, 0, static_cast<int>(map.width) - 1);
            cy = std::clamp(cy, 0, static_cast<int>(map.height) - 1);
            return layers[static_cast<usize>(cy) * map.width + cx];
        };
        const auto sample_height = [&](float wx, float wy) -> float {
            if (wx < 0.0F || wy < 0.0F || wx >= world_w || wy >= world_h) return -1.0e9F;
            const auto cx = std::min(map.width - 1U, static_cast<uint32>(wx / cell_size));
            const auto cy = std::min(map.height - 1U, static_cast<uint32>((world_h - wy) / cell_size));
            return height[static_cast<usize>(cy) * map.width + cx];
        };

        const auto mix_color = [](uint32 a, uint32 b, float w) -> uint32 { return detail::mix_color(a, b, w); };

        float lx = options.sun_dir[0];
        float ly = options.sun_dir[1];
        float lz = options.sun_dir[2];
        const auto llen = std::sqrt(lx * lx + ly * ly + lz * lz);
        lx /= llen;
        ly /= llen;
        lz /= llen;

        const auto target_z = sample_height(camera.target_x, camera.target_y);
        const auto ground = target_z > -1.0e8F ? target_z : 0.0F;
        const auto pitch = std::clamp(camera.pitch, 0.15F, 1.45F);
        const auto dist = camera.height / std::sin(pitch);
        const float fx = std::cos(pitch) * std::sin(camera.yaw);
        const float fy = std::cos(pitch) * std::cos(camera.yaw);
        const float fz = -std::sin(pitch);
        auto rx = fy;
        auto ry = -fx;
        auto rz = 0.0F;
        const auto rlen = std::sqrt(rx * rx + ry * ry + rz * rz);
        rx /= rlen;
        ry /= rlen;
        rz /= rlen;
        auto ux = ry * fz - rz * fy;
        auto uy = rz * fx - rx * fz;
        auto uz = rx * fy - ry * fx;
        const auto ulen = std::sqrt(ux * ux + uy * uy + uz * uz);
        ux /= ulen;
        uy /= ulen;
        uz /= ulen;
        const auto cam_x = camera.target_x - fx * dist;
        const auto cam_y = camera.target_y - fy * dist;
        const auto cam_z = ground + camera.height;
        const auto tan_half = std::tan(camera.fov * 0.5F);
        const auto aspect = static_cast<float>(rw) / static_cast<float>(rh);
        const auto sun_horiz = std::max(1.0e-4F, std::sqrt(lx * lx + ly * ly));
        const auto step_x = lx / sun_horiz;
        const auto step_y = ly / sun_horiz;

        const auto hit_color = [&](float hx, float hy, float hz) -> uint32 {
            if (map.has_water && hz <= map.water_plane_z + 0.01F) return options.water_color;
            const auto cx = std::clamp(static_cast<int>(hx / cell_size), 0, static_cast<int>(map.width) - 1);
            const auto cy = std::clamp(static_cast<int>((world_h - hy) / cell_size), 0, static_cast<int>(map.height) - 1);
            const auto fxc = hx / cell_size - static_cast<float>(cx);
            const auto fyc = (world_h - hy) / cell_size - static_cast<float>(cy);
            const auto wx = hx / cell_size;
            const auto wy = (world_h - hy) / cell_size;
            const auto &layers = cell_at(cx, cy);
            const auto dir1 = static_cast<uint32>(layers.packed & 0xFU);
            const auto flags1 = static_cast<uint32>((layers.packed >> 4U) & 0x3U);
            const auto dir2 = static_cast<uint32>((layers.packed >> 8U) & 0xFU);
            const auto flags2 = static_cast<uint32>((layers.packed >> 12U) & 0x3U);
            const auto blend1 = detail::blend_factor(dir1, flags1, fxc, fyc);
            const auto blend2 = detail::blend_factor(dir2, flags2, fxc, fyc);
            const auto c0 = detail::sample_layer_bilinear(layers.base, wx, wy);
            const auto c1 = blend1 > 0.0F ? detail::sample_layer_bilinear(layers.blend, wx, wy) : c0;
            const auto c2 = blend2 > 0.0F ? detail::sample_layer_bilinear(layers.three, wx, wy) : c0;
            auto base = mix_color(mix_color(c0, c1, blend1), c2, blend2);

            const auto h_l = sample_height(hx - cell_size, hy);
            const auto h_r = sample_height(hx + cell_size, hy);
            const auto h_d = sample_height(hx, hy - cell_size);
            const auto h_u = sample_height(hx, hy + cell_size);
            const auto nx = h_l - h_r;
            const auto ny = h_d - h_u;
            const auto nz = 2.0F * cell_size;
            const auto nlen = std::sqrt(nx * nx + ny * ny + nz * nz);
            auto lambert = std::max(0.0F, (nx / nlen) * lx + (ny / nlen) * ly + (nz / nlen) * lz);
            for (uint32 s = 1; s <= options.shadow_steps; ++s) {
                const auto sx = hx + step_x * static_cast<float>(s) * cell_size;
                const auto sy = hy + step_y * static_cast<float>(s) * cell_size;
                const auto ray_z = hz + (lz / sun_horiz) * static_cast<float>(s) * cell_size;
                if (sample_height(sx, sy) > ray_z) {
                    lambert *= options.shadow_strength;
                    break;
                }
            }
            const auto shade = std::clamp(options.ambient + (1.0F - options.ambient) * lambert, 0.0F, 1.4F);
            const auto mod = [&](uint32 c) { return static_cast<uint32>(std::clamp(static_cast<float>(c) * shade, 0.0F, 255.0F)); };
            return argb(static_cast<uint8>(mod((base >> 16U) & 0xFFU)), static_cast<uint8>(mod((base >> 8U) & 0xFFU)), static_cast<uint8>(mod(base & 0xFFU)));
        };

        image hi(rw, rh, argb(0, 0, 0));
        std::vector<float> zbuf(static_cast<usize>(rw) * rh, 1.0e30F);
        for (uint32 py = 0; py < rh; ++py) {
            const auto ndc_y = 1.0F - 2.0F * (static_cast<float>(py) + 0.5F) / static_cast<float>(rh);
            for (uint32 px = 0; px < rw; ++px) {
                const auto ndc_x = 2.0F * (static_cast<float>(px) + 0.5F) / static_cast<float>(rw) - 1.0F;
                float dx = fx + rx * ndc_x * tan_half * aspect + ux * ndc_y * tan_half;
                float dy = fy + ry * ndc_x * tan_half * aspect + uy * ndc_y * tan_half;
                float dz = fz + rz * ndc_x * tan_half * aspect + uz * ndc_y * tan_half;
                const auto dlen = std::sqrt(dx * dx + dy * dy + dz * dz);
                dx /= dlen;
                dy /= dlen;
                dz /= dlen;
                const auto pixel = static_cast<usize>(py) * rw + px;

                if (dz >= -1.0e-4F) {
                    const auto t = std::clamp(std::abs(dz) * 2.0F, 0.0F, 1.0F);
                    hi.data()[pixel] =
                        argb(static_cast<uint8>(120.0F + 60.0F * t), static_cast<uint8>(150.0F + 45.0F * t), static_cast<uint8>(190.0F + 40.0F * t));
                    continue;
                }

                auto t = cell_size * 0.5F;
                auto dt = cell_size * 0.5F;
                auto prev_t = t;
                bool hit = false;
                float hit_t = 0.0F;
                for (int iter = 0; iter < 4000 && t < 60000.0F; ++iter) {
                    const auto wx = cam_x + dx * t;
                    const auto wy = cam_y + dy * t;
                    const auto wz = cam_z + dz * t;
                    if (wx < 0.0F || wy < 0.0F || wx >= world_w || wy >= world_h) {
                        prev_t = t;
                        dt *= 1.03F;
                        t += dt;
                        continue;
                    }
                    if (map.has_water && wz <= map.water_plane_z) {
                        hit = true;
                        hit_t = t;
                        break;
                    }
                    if (wz <= sample_height(wx, wy)) {
                        hit = true;
                        hit_t = t;
                        break;
                    }
                    prev_t = t;
                    dt *= 1.03F;
                    t += dt;
                }
                if (!hit) {
                    hi.data()[pixel] = argb(150, 170, 200);
                    continue;
                }

                auto lo = prev_t;
                auto up = hit_t;
                for (int i = 0; i < 6; ++i) {
                    const auto mid = 0.5F * (lo + up);
                    const auto wx = cam_x + dx * mid;
                    const auto wy = cam_y + dy * mid;
                    const auto wz = cam_z + dz * mid;
                    const auto water = map.has_water && wz <= map.water_plane_z;
                    if (water || wz <= sample_height(wx, wy)) {
                        up = mid;
                    } else {
                        lo = mid;
                    }
                }
                hi.data()[pixel] = hit_color(cam_x + dx * up, cam_y + dy * up, cam_z + dz * up);
                // View-space depth of the hit, for the object pass below.
                zbuf[pixel] = up * (dx * fx + dy * fy + dz * fz);
            }
        }

        if (objects != nullptr && !objects->empty()) {
            detail::rasterize_objects(hi, zbuf, *objects, {cam_x, cam_y, cam_z}, {fx, fy, fz}, {rx, ry, rz}, {ux, uy, uz}, tan_half, aspect);
        }

        if (ss == 1U) return hi;
        image out(out_w, out_h, argb(0, 0, 0));
        for (uint32 y = 0; y < out_h; ++y) {
            for (uint32 x = 0; x < out_w; ++x) {
                uint32 r = 0;
                uint32 g = 0;
                uint32 b = 0;
                for (uint32 sy = 0; sy < ss; ++sy) {
                    for (uint32 sx = 0; sx < ss; ++sx) {
                        const auto c = hi.data()[static_cast<usize>(y * ss + sy) * rw + x * ss + sx];
                        r += (c >> 16U) & 0xFFU;
                        g += (c >> 8U) & 0xFFU;
                        b += c & 0xFFU;
                    }
                }
                const auto n = ss * ss;
                out.data()[static_cast<usize>(y) * out_w + x] = argb(static_cast<uint8>(r / n), static_cast<uint8>(g / n), static_cast<uint8>(b / n));
            }
        }
        return out;
    }
}