"""Pyramide de tuiles XYZ (EPSG:3857) rendue à la demande depuis les dalles.

Schéma de tuilage identique à celui d'OpenStreetMap / Google Maps : grille
Web Mercator, origine au coin nord-ouest, `{z}/{x}/{y}`, tuiles de 256 px
(512 px avec `scale=2`, convention `@2x`). Les rendus du pipeline restent des
dalles Lambert 93 de 1 km : chaque tuile est composée à la volée en reprojetant
les dalles qui l'intersectent (transformation projective par dalle, erreur très
inférieure au pixel), puis mise en cache sur disque.

Volontairement sans GDAL ni numpy : Pillow + pyproj suffisent, l'image légère
(Dockerfile.maps / Dockerfile.webapp) reste petite et portable ARM64.
"""

import logging
import math
import os
import re
import threading
import time
from collections import OrderedDict
from functools import lru_cache
from pathlib import Path

logger = logging.getLogger("lidar")

# --- Contrat de tuilage (cf. docs/MAPS.md) ---------------------------------
TILE_SIZE = 256                # taille canonique (OSM/XYZ) ; @2x → 512
TILE_MIN_Z = 5
TILE_MAX_NATIVE_Z = 19         # 0,2 m/px ≈ résolution du z19 à la latitude 47°
TILE_DIRNAME = "index_xyz"     # cache disque, sous le dossier de sortie
WEBP_QUALITY = 78
AVIF_QUALITY = 60
# PNG palettisé (PNG8 + alpha) : ~5× plus léger (170 → 32 Ko sur une dalle
# réelle) pour un écart moyen de ~4 niveaux sur une rampe de couleur. Laissé
# DÉSACTIVÉ par défaut : le PNG canonique reste sans perte, la fidélité prime
# sur le débit pour un produit d'interprétation. `LIDAR_TILE_PNG_PALETTE=1`
# l'active quand la bande passante compte (consultation mobile).
PNG_PALETTE = os.environ.get("LIDAR_TILE_PNG_PALETTE", "") == "1"

# Demi-circonférence équatoriale : emprise du Web Mercator (EPSG:3857).
ORIGIN = 20037508.342789244

# Résolutions nominales des paliers de source réutilisés tels quels par le
# pipeline (m/px) : vignette 256 px/km, vignette intermédiaire 640 px/km.
_THUMB_RES = 1000.0 / 256
_MID_RES = 1000.0 / 640

_SUBTILE_RE = re.compile(r"_(\d+)_(\d+)\.avif$")

# Serveur de dalles amont (webapp du pipeline) : quand il est défini, l'index
# des sources vient de son /api/tiles et les images manquantes sont rapatriées
# à la demande dans le cache local — le conteneur carte n'a alors besoin
# d'aucune donnée locale au démarrage.
REMOTE_SOURCE_URL = (os.environ.get("LIDAR_SOURCE_URL") or "").rstrip("/")
REMOTE_SOURCE_TOKEN = os.environ.get("LIDAR_SOURCE_TOKEN") or None
_REMOTE_TTL = 60.0
_remote_cache = {"payload": None, "at": 0.0, "index": None, "root": None}
# Téléchargements de dalles simultanés (réseau) — modeste par défaut : sur un
# petit serveur, chaque source rapatriée est un fichier de plusieurs Mo.
FETCH_WORKERS = max(1, int(os.environ.get("LIDAR_TILE_FETCH_WORKERS", "2") or 2))
_fetch_sem = threading.Semaphore(FETCH_WORKERS)
_fetch_locks = {}
_fetch_guard = threading.Lock()

# Index des sources reconstruit au plus toutes les _INDEX_TTL secondes (ou dès
# qu'un dossier change de mtime : ajout/suppression de dalle).
_INDEX_TTL = 20.0
_index_cache = {}


# ---------------------------------------------------------------------------
# Géométrie de la grille
# ---------------------------------------------------------------------------

def tile_bounds_3857(z, x, y):
    """Emprise (ouest, sud, est, nord) d'une tuile XYZ en mètres EPSG:3857."""
    span = 2.0 * ORIGIN / (2 ** z)
    west = -ORIGIN + x * span
    north = ORIGIN - y * span
    return west, north - span, west + span, north


def ground_resolution(z, scale=1):
    """Résolution d'une tuile à l'équateur (m/px) ; ×cos(lat) sur le terrain."""
    return 2.0 * ORIGIN / (TILE_SIZE * scale * (2 ** z))


def tile_latitude(z, y):
    """Latitude (degrés) du centre d'une tuile — sert au calage de résolution."""
    n = math.pi - 2.0 * math.pi * (y + 0.5) / (2 ** z)
    return math.degrees(math.atan(math.sinh(n)))


def target_resolution(z, y, scale=1):
    """Résolution terrain visée par la tuile (m/px), latitude comprise."""
    return ground_resolution(z, scale) * math.cos(math.radians(tile_latitude(z, y)))


@lru_cache(maxsize=4)
def _transformer(src, dst):
    from pyproj import Transformer
    return Transformer.from_crs(src, dst, always_xy=True)


def to_l93(xs, ys):
    """EPSG:3857 → EPSG:2154 (listes de coordonnées)."""
    return _transformer("EPSG:3857", "EPSG:2154").transform(xs, ys)


def to_3857(xs, ys):
    """EPSG:2154 → EPSG:3857 (listes de coordonnées)."""
    return _transformer("EPSG:2154", "EPSG:3857").transform(xs, ys)


@lru_cache(maxsize=4096)
def wgs84_to_l93(lon, lat):
    """Point WGS84 → Lambert 93 (utilisé par la fiche d'information dalle)."""
    return _transformer("EPSG:4326", "EPSG:2154").transform(lon, lat)


@lru_cache(maxsize=4096)
def tile_bounds_l93(z, x, y, samples=5):
    """Emprise L93 englobant une tuile XYZ.

    Les bords d'une tuile ne sont pas des droites en Lambert 93 : on
    échantillonne une grille `samples`×`samples` plutôt que les seuls coins,
    sinon l'emprise est sous-estimée aux petits zooms (tuiles de centaines de km).
    """
    west, south, east, north = tile_bounds_3857(z, x, y)
    xs, ys = [], []
    for i in range(samples):
        for j in range(samples):
            xs.append(west + (east - west) * i / (samples - 1))
            ys.append(south + (north - south) * j / (samples - 1))
    lx, ly = to_l93(xs, ys)
    return min(lx), min(ly), max(lx), max(ly)


# ---------------------------------------------------------------------------
# Transformation projective (mapping sortie → source, convention Pillow)
# ---------------------------------------------------------------------------

def _solve(matrix, rhs):
    """Résout un système linéaire dense (pivot partiel), sans numpy."""
    n = len(rhs)
    a = [row[:] + [rhs[i]] for i, row in enumerate(matrix)]
    for col in range(n):
        piv = max(range(col, n), key=lambda r: abs(a[r][col]))
        if abs(a[piv][col]) < 1e-12:
            return None
        a[col], a[piv] = a[piv], a[col]
        inv = 1.0 / a[col][col]
        for k in range(col, n + 1):
            a[col][k] *= inv
        for r in range(n):
            if r == col:
                continue
            f = a[r][col]
            if f:
                for k in range(col, n + 1):
                    a[r][k] -= f * a[col][k]
    return [a[i][n] for i in range(n)]


def perspective_coeffs(dst_quad, src_quad):
    """Coefficients Pillow `Image.PERSPECTIVE` mappant sortie → source.

    Pillow échantillonne la SOURCE en (x', y') = ((a x + b y + c) / (g x + h y + 1),
    (d x + e y + f) / (g x + h y + 1)) pour chaque pixel (x, y) de la SORTIE :
    on résout donc les 8 inconnues à partir de 4 correspondances
    (point de sortie → point source). Retourne None si le quadrilatère est
    dégénéré (dalle réduite à un point au dézoom extrême).
    """
    matrix, rhs = [], []
    for (dx, dy), (sx, sy) in zip(dst_quad, src_quad):
        matrix.append([dx, dy, 1, 0, 0, 0, -sx * dx, -sx * dy])
        rhs.append(sx)
        matrix.append([0, 0, 0, dx, dy, 1, -sy * dx, -sy * dy])
        rhs.append(sy)
    return _solve(matrix, rhs)


# ---------------------------------------------------------------------------
# Index des sources : dalles du pipeline, par couche et par palier
# ---------------------------------------------------------------------------

class _Source:
    """Une image source géoréférencée (dalle entière ou quadrant).

    `url` non nul : l'image vit sur le serveur de dalles amont et n'est
    rapatriée qu'au premier besoin (`ensure()`), dans `path`. `version` est
    alors la mtime (ms) annoncée par l'amont : elle sert de date de référence
    pour la péremption des tuiles, avant même tout téléchargement.
    """

    __slots__ = ("path", "bounds", "res", "url", "version")

    def __init__(self, path, bounds, res, url=None, version=None):
        self.path = path
        self.bounds = bounds      # (min_x, min_y, max_x, max_y) en L93
        self.res = res            # résolution nominale (m/px)
        self.url = url
        self.version = version

    def mtime(self):
        if self.version is not None:
            return self.version / 1000.0
        try:
            return self.path.stat().st_mtime
        except OSError:
            return None

    def ensure(self):
        """Garantit la présence locale de l'image (rapatriement si besoin)."""
        if self.url is None or self.path.is_file():
            return self.path.is_file()
        return _fetch_source(self.url, self.path)


def _remote_thumb_px(k):
    """Côté (px) de la vignette servie par l'amont : 256 par dalle, 160 par quadrant."""
    return 160.0 if k > 1 else 256.0


def _fetch_source(url, dest):
    """Télécharge une image source depuis l'amont (écriture atomique)."""
    import urllib.request
    with _fetch_guard:
        lock = _fetch_locks.setdefault(str(dest), threading.Lock())
    with lock:
        if dest.is_file():
            return True
        try:
            req = urllib.request.Request(
                url, headers={"User-Agent": "lidar-maps-source"})
            if REMOTE_SOURCE_TOKEN:
                req.add_header("X-Lidar-Token", REMOTE_SOURCE_TOKEN)
            with _fetch_sem, urllib.request.urlopen(req, timeout=120) as r:
                data = r.read()
        except Exception as e:  # noqa: BLE001 — amont éteint : tuile partielle
            logger.debug(f"Source amont indisponible ({url}) : {e}")
            return False
        _write_atomic(dest, data)
        logger.info(f"Source rapatriée : {dest.name} ({len(data) / 1e6:.1f} Mo)")
        return dest.is_file()


def _remote_payload(force=False):
    """Index du serveur de dalles amont (/api/tiles), en cache 60 s."""
    import json
    import urllib.request
    now = time.time()
    if not force and _remote_cache["payload"] is not None \
            and now - _remote_cache["at"] < _REMOTE_TTL:
        return _remote_cache["payload"]
    try:
        req = urllib.request.Request(f"{REMOTE_SOURCE_URL}/api/tiles",
                                     headers={"User-Agent": "lidar-maps-source"})
        with urllib.request.urlopen(req, timeout=30) as r:
            payload = json.loads(r.read().decode("utf-8"))
    except Exception as e:  # noqa: BLE001 — on garde le dernier index connu
        logger.warning(f"Index amont injoignable ({REMOTE_SOURCE_URL}) : {e}")
        payload = _remote_cache["payload"]
    _remote_cache.update(payload=payload, at=now)
    return payload


def _remote_index(output_dir, force=False):
    """Inventaire construit depuis l'amont : sources non encore rapatriées.

    Le résultat est mémoïsé avec la charge utile : le reconstruire coûte des
    centaines de millisecondes sur un catalogue de plusieurs milliers de
    dalles, ce qui se paierait à CHAQUE tuile servie.
    """
    payload = _remote_payload(force)
    if not payload or not payload.get("tiles"):
        return {}
    if (_remote_cache["index"] is not None
            and _remote_cache["root"] == str(output_dir)
            and _remote_cache.get("built") is payload):
        return _remote_cache["index"]
    output_dir = Path(output_dir)
    layers = {}
    for entry in payload["tiles"]:
        col, row = entry.get("col"), entry.get("row")
        if col is None or row is None:
            continue
        res = float(entry.get("resolution") or 0.5)
        k = int(entry.get("sub_k") or 1)
        step = 1000.0 / k
        i, j = int(entry.get("sub_i") or 0), int(entry.get("sub_j") or 0)
        base = _cell_bounds(col, row)
        bounds = (base[0] + i * step, base[1] + j * step,
                  base[0] + (i + 1) * step, base[1] + (j + 1) * step)
        for viz, info in (entry.get("viz") or {}).items():
            by_res = layers.setdefault(viz, {}).setdefault((col, row), {})
            for tier, nominal in (("thumb", step / _remote_thumb_px(k)),
                                  ("mid", step / 640.0),
                                  ("full", res)):
                url_rel = info.get(tier)
                if not url_rel:
                    continue
                rel, _, query = url_rel.partition("?")
                version = None
                if query.startswith("v="):
                    try:
                        version = int(query[2:])
                    except ValueError:
                        version = None
                rank = 0 if k > 1 else 1
                by_res.setdefault((nominal, rank), []).append(_Source(
                    output_dir / rel, bounds, nominal,
                    url=f"{REMOTE_SOURCE_URL}/{url_rel}", version=version))
    out = {}
    for viz, per_cell in layers.items():
        for cell, by_res in per_cell.items():
            out.setdefault(viz, {})[cell] = [
                by_res[key] for key in sorted(by_res, key=lambda t: (-t[0], t[1]))]
    _remote_cache.update(index=out, root=str(output_dir), built=payload)
    return out


def _cell_bounds(col, row):
    """Emprise L93 d'une dalle LHD : X ∈ [col, col+1] km, Y ∈ [row-1, row] km."""
    return (col * 1000.0, (row - 1) * 1000.0, (col + 1) * 1000.0, row * 1000.0)


def cell_bounds_wgs84(col, row):
    """Emprise WGS84 [ouest, sud, est, nord] d'une dalle LHD (4 coins L93)."""
    min_x, min_y, max_x, max_y = _cell_bounds(col, row)
    xs = (min_x, min_x, max_x, max_x)
    ys = (min_y, max_y, min_y, max_y)
    lons, lats = _transformer("EPSG:2154", "EPSG:4326").transform(xs, ys)
    return (min(lons), min(lats), max(lons), max(lats))


def _dir_mtime(path):
    try:
        return path.stat().st_mtime
    except OSError:
        return None


def _split_viz_key(stem, known):
    """`{dir_name}_{viz}` → (dir_name, viz) en s'appuyant sur les clés connues.

    Les clés de visualisation contiennent des soulignés (`positive_openness`,
    `hillshade_multi`) : on ne peut pas couper au dernier `_`, on reconnaît un
    suffixe connu (le plus long d'abord).
    """
    for key in known:
        if stem.endswith("_" + key):
            return stem[:-len(key) - 1], key
    return None, None


def _cell_of_dir(dir_name):
    """(col, row, résolution) d'un nom de dossier de dalle, ou None."""
    from .index import _strip_res_suffix, parse_basename_coords
    coords = parse_basename_coords(dir_name)
    if coords is None:
        return None
    _base, res = _strip_res_suffix(dir_name)
    return coords[0], coords[1], res


def _build_index(output_dir):
    """Inventaire {couche: {(col, row): [paliers du plus grossier au plus fin]}}.

    Les trois sources du pipeline sont scannées INDÉPENDAMMENT — dalles
    (`visualisations/`), quadrants (`index_subtiles/`) et vignettes
    (`index_thumbs/`). Un cache partiel reste donc exploitable : sur une
    machine légère, seuls quadrants et vignettes sont rapatriés, jamais les
    dalles entières.
    """
    # _SUBTILE_THUMB_PX : taille des vignettes de quadrant, définie par
    # l'index (la changer là-bas doit rester sans effet ici).
    from .index import _SUBTILE_THUMB_PX, VIZ_LABELS, scan_tiles

    output_dir = Path(output_dir)
    vis_dir = output_dir / "visualisations"
    thumb_dir = output_dir / "index_thumbs"
    sub_dir = output_dir / "index_subtiles"

    # res → sources, par couche et par dalle ; trié en paliers à la fin.
    records = {}

    def add(layer, col, row, res, source, rank=1):
        # Clé de palier = (résolution, rang) : à résolution égale, les
        # quadrants (rang 0) passent avant la dalle entière (rang 1) — même
        # rendu, 4× moins de pixels à décoder.
        records.setdefault(layer, {}).setdefault((col, row), {}) \
            .setdefault((res, rank), []).append(source)

    # 1. Dalles entières (palier le plus fin quand il est présent).
    known = set(VIZ_LABELS)
    for tile in scan_tiles(vis_dir):
        col, row, res = tile["col"], tile["row"], tile["resolution"]
        cell = _cell_bounds(col, row)
        dir_path = Path(tile["dir_path"])
        for viz_key, info in tile["viz"].items():
            known.add(viz_key)
            source_dir = dir_path
            if info.get("dir_name") and info["dir_name"] != dir_path.name:
                source_dir = dir_path.parent / info["dir_name"]
            full = source_dir / info["filename"]
            if full.is_file():
                add(viz_key, col, row, res, _Source(full, cell, res))

    # Clés les plus longues d'abord : `positive_openness` avant `openness`.
    known = sorted(known, key=len, reverse=True)

    # 2. Quadrants (index_subtiles) : AVIF pleine résolution + ses vignettes.
    if sub_dir.is_dir():
        quads = {}
        for f in sub_dir.iterdir():
            name = f.name
            for suffix, tier in ((".avif", "full"), ("_mid.webp", "mid"),
                                 (f"_thumb{_SUBTILE_THUMB_PX}.webp", "thumb")):
                if not name.endswith(suffix):
                    continue
                m = _SUBTILE_RE.search(name[:-len(suffix)] + ".avif")
                if not m:
                    break
                stem = name[:m.start()]
                i, j = int(m.group(1)), int(m.group(2))
                dir_name, viz = _split_viz_key(stem, known)
                if dir_name is None:
                    break
                quads.setdefault((dir_name, viz, tier), []).append((i, j, f))
                break
        for (dir_name, viz, tier), items in quads.items():
            cell = _cell_of_dir(dir_name)
            if cell is None:
                continue
            col, row, res = cell
            k = max(max(i for i, _j, _f in items),
                    max(j for _i, j, _f in items)) + 1
            step = 1000.0 / k
            px = {"full": step / res, "mid": 640, "thumb": _SUBTILE_THUMB_PX}[tier]
            tier_res = res if tier == "full" else step / px
            base = _cell_bounds(col, row)
            for i, j, f in items:
                add(viz, col, row, tier_res, _Source(
                    f, (base[0] + i * step, base[1] + j * step,
                        base[0] + (i + 1) * step, base[1] + (j + 1) * step),
                    tier_res), rank=0)

    # 3. Vignettes de dalle (paliers grossiers).
    if thumb_dir.is_dir():
        for f in thumb_dir.iterdir():
            if f.suffix.lower() not in (".jpg", ".jpeg", ".webp", ".png"):
                continue
            stem = f.name[:-len(f.suffix)]
            res_px = _MID_RES if stem.endswith("_mid") else _THUMB_RES
            if stem.endswith("_mid"):
                stem = stem[:-4]
            dir_name, viz = _split_viz_key(stem, known)
            if dir_name is None:
                continue
            cell = _cell_of_dir(dir_name)
            if cell is None:
                continue
            col, row, _res = cell
            add(viz, col, row, res_px, _Source(f, _cell_bounds(col, row), res_px))

    # Paliers triés du plus grossier au plus fin (l'ordre de choix du rendu).
    layers = {}
    for viz, per_cell in records.items():
        for (col, row), by_res in per_cell.items():
            tiers = [by_res[k] for k in sorted(by_res, key=lambda k: (-k[0], k[1]))]
            layers.setdefault(viz, {})[(col, row)] = tiers
    return layers


def source_index(output_dir, force=False):
    """Index des sources, mémoïsé (TTL + mtime des dossiers surveillés)."""
    output_dir = Path(output_dir)
    key = str(output_dir)
    stamp = tuple(_dir_mtime(output_dir / d) for d in
                  ("visualisations", "index_thumbs", "index_subtiles"))
    entry = _index_cache.get(key)
    now = time.time()
    if entry and not force and now - entry["at"] < _INDEX_TTL:
        # Le TTL prime sur la mtime des dossiers : en mode amont, chaque source
        # rapatriée la modifierait et provoquerait un rescan par tuile servie.
        if REMOTE_SOURCE_URL or entry["stamp"] == stamp:
            return entry["layers"]
    if REMOTE_SOURCE_URL:
        # L'amont fait autorité : il connaît toutes les dalles, le cache local
        # n'en détient qu'une partie (et grossit à chaque source rapatriée —
        # le rescanner à chaque tuile coûterait plus cher que le rendu). Les
        # sources déjà présentes sont servies depuis le disque (_Source.ensure).
        layers = dict(_remote_index(output_dir, force))
        if not layers:
            layers = _build_index(output_dir)   # amont muet : cache local seul
    else:
        layers = _build_index(output_dir)
    _index_cache[key] = {"layers": layers, "stamp": stamp, "at": now}
    return layers


def available_layers(output_dir):
    """Couches présentes sur disque, ordonnées comme le panneau de la carte."""
    from .index import _VIZ_FALLBACK_ORDER
    found = set(source_index(output_dir))
    ordered = [v for v in _VIZ_FALLBACK_ORDER if v in found]
    ordered += sorted(found - set(ordered))
    return ordered


def grid_bounds_l93(output_dir, layer=None):
    """Emprise L93 (min_x, min_y, max_x, max_y) des dalles disponibles."""
    layers = source_index(output_dir)
    cells = set()
    for key, per_cell in layers.items():
        if layer and key != layer:
            continue
        cells.update(per_cell)
    if not cells:
        return None
    xs = [c for c, _r in cells]
    ys = [r for _c, r in cells]
    return (min(xs) * 1000.0, (min(ys) - 1) * 1000.0,
            (max(xs) + 1) * 1000.0, max(ys) * 1000.0)


def grid_bounds_wgs84(output_dir, layer=None):
    """Emprise WGS84 [ouest, sud, est, nord] des dalles (TileJSON, WMTS)."""
    l93 = grid_bounds_l93(output_dir, layer)
    if l93 is None:
        return None
    min_x, min_y, max_x, max_y = l93
    tr = _transformer("EPSG:2154", "EPSG:4326")
    xs = [min_x, max_x, min_x, max_x]
    ys = [min_y, min_y, max_y, max_y]
    lons, lats = tr.transform(xs, ys)
    return [min(lons), min(lats), max(lons), max(lats)]


def tiles_stamp(output_dir):
    """Version globale du jeu de tuiles (max des mtimes) pour l'URL du client."""
    newest = 0.0
    for per_cell in source_index(output_dir).values():
        for tiers in per_cell.values():
            for src in tiers[-1]:
                m = src.mtime()
                if m and m > newest:
                    newest = m
    return int(newest * 1000)


# ---------------------------------------------------------------------------
# Rendu d'une tuile
# ---------------------------------------------------------------------------

# Cache des images sources décodées : le décodage (AVIF surtout) domine le
# coût d'une tuile et les tuiles voisines partagent leurs dalles. Le budget est
# exprimé en OCTETS, pas en nombre d'entrées : une dalle 5000² pèse ~75 Mo
# quand une vignette en pèse 0,2 — un cache « N entrées » ferait déborder la
# mémoire d'une petite machine.
SOURCE_CACHE_BYTES = int(os.environ.get("LIDAR_TILE_SOURCE_CACHE_MB", "192")) * 1024 * 1024
_source_cache = OrderedDict()
_source_cache_lock = threading.Lock()


def _open_source(path_str, mtime):
    """Image source décodée, mémoïsée par (chemin, mtime).

    `mtime` fait partie de la clé : une dalle régénérée invalide l'entrée.
    """
    from PIL import Image
    key = (str(path_str), mtime)
    with _source_cache_lock:
        img = _source_cache.get(key)
        if img is not None:
            _source_cache.move_to_end(key)
            return img
    img = Image.open(str(path_str))
    img.load()
    if img.mode not in ("RGB", "RGBA"):
        img = img.convert("RGB")
    size = img.size[0] * img.size[1] * len(img.getbands())
    with _source_cache_lock:
        _source_cache[key] = img
        _source_cache[key].info["_bytes"] = size
        total = sum(i.info.get("_bytes", 0) for i in _source_cache.values())
        while total > SOURCE_CACHE_BYTES and len(_source_cache) > 1:
            _k, old_img = _source_cache.popitem(last=False)
            total -= old_img.info.get("_bytes", 0)
    return img


def clear_source_cache():
    """Vide le cache d'images sources (tests, pression mémoire)."""
    with _source_cache_lock:
        _source_cache.clear()


def _pick_tier(tiers, target_res):
    """Palier le plus grossier dont la résolution suffit à la tuile visée."""
    for group in tiers:
        if group[0].res <= target_res:
            return group
    return tiers[-1]


def _contributing(output_dir, layer, z, x, y, scale):
    """Sources intersectant la tuile, palier choisi selon la résolution visée."""
    per_cell = source_index(output_dir).get(layer)
    if not per_cell:
        return []
    min_x, min_y, max_x, max_y = tile_bounds_l93(z, x, y)
    target = target_resolution(z, y, scale)
    out = []
    for (col, row), tiers in per_cell.items():
        b = _cell_bounds(col, row)
        if b[2] <= min_x or b[0] >= max_x or b[3] <= min_y or b[1] >= max_y:
            continue
        for src in _pick_tier(tiers, target):
            s = src.bounds
            if s[2] <= min_x or s[0] >= max_x or s[3] <= min_y or s[1] >= max_y:
                continue
            out.append(src)
    return out


def _paste_source(canvas, src, z, x, y, size, resample):
    """Reprojette une source dans la tuile (transformation projective)."""
    from PIL import Image

    if not src.ensure():
        return False
    mtime = src.mtime()
    if mtime is None:
        return False
    try:
        img = _open_source(str(src.path), mtime)
    except Exception as e:  # noqa: BLE001 — source illisible : tuile partielle
        logger.debug(f"Source de tuile illisible ({src.path.name}) : {e}")
        return False

    w, h = img.size
    min_x, min_y, max_x, max_y = src.bounds
    res_x = (max_x - min_x) / w
    res_y = (max_y - min_y) / h

    # Fenêtre source utile = intersection avec l'emprise L93 de la tuile,
    # élargie de 2 px pour que l'interpolation dispose de son voisinage.
    t_min_x, t_min_y, t_max_x, t_max_y = tile_bounds_l93(z, x, y)
    c0 = max(0, int(math.floor((max(t_min_x, min_x) - min_x) / res_x)) - 2)
    c1 = min(w, int(math.ceil((min(t_max_x, max_x) - min_x) / res_x)) + 2)
    r0 = max(0, int(math.floor((max_y - min(t_max_y, max_y)) / res_y)) - 2)
    r1 = min(h, int(math.ceil((max_y - max(t_min_y, min_y)) / res_y)) + 2)
    if c1 - c0 < 1 or r1 - r0 < 1:
        return False

    crop = img.crop((c0, r0, c1, r1))
    if crop.mode != "RGBA":
        crop = crop.convert("RGBA")
    cw, ch = crop.size

    # Coins L93 de la fenêtre découpée → EPSG:3857 → pixels de la tuile.
    wx0 = min_x + c0 * res_x
    wx1 = min_x + c1 * res_x
    wy1 = max_y - r0 * res_y
    wy0 = max_y - r1 * res_y
    mx, my = to_3857([wx0, wx1, wx1, wx0], [wy1, wy1, wy0, wy0])
    west, south, east, north = tile_bounds_3857(z, x, y)
    px = [(m - west) / (east - west) * size for m in mx]
    py = [(north - m) / (north - south) * size for m in my]
    dst_quad = list(zip(px, py))                       # NW, NE, SE, SW
    src_quad = [(0, 0), (cw, 0), (cw, ch), (0, ch)]
    coeffs = perspective_coeffs(dst_quad, src_quad)
    if coeffs is None:
        return False

    warped = crop.transform((size, size), Image.PERSPECTIVE, coeffs,
                            resample=resample, fillcolor=(0, 0, 0, 0))
    canvas.alpha_composite(warped)
    return True


def render_tile(output_dir, layer, z, x, y, scale=1):
    """Rend une tuile en mémoire. Retourne une image RGBA, ou None si vide."""
    from PIL import Image

    sources = _contributing(output_dir, layer, z, x, y, scale)
    if not sources:
        return None
    size = TILE_SIZE * scale
    canvas = Image.new("RGBA", (size, size), (0, 0, 0, 0))
    # Agrandissement (zoom natif) : bicubique ; réduction : bilinéaire suffit
    # puisque le palier source est déjà calé sur la résolution de la tuile.
    resample = Image.BICUBIC
    painted = False
    for src in sources:
        painted |= _paste_source(canvas, src, z, x, y, size, resample)
    if not painted:
        return None
    return canvas


# ---------------------------------------------------------------------------
# Cache disque
# ---------------------------------------------------------------------------

def tile_cache_path(output_dir, layer, z, x, y, scale=1, fmt="png"):
    """Chemin du fichier de cache d'une tuile."""
    suffix = "@2x" if scale == 2 else ""
    return (Path(output_dir) / TILE_DIRNAME / layer / str(z) / str(x)
            / f"{y}{suffix}.{fmt}")


def _encode(img, fmt):
    import io
    buf = io.BytesIO()
    if fmt == "png":
        if PNG_PALETTE:
            from PIL import Image
            img = img.quantize(colors=255, method=Image.Quantize.FASTOCTREE)
        img.save(buf, format="PNG", optimize=False)
    elif fmt == "webp":
        img.save(buf, format="WEBP", quality=WEBP_QUALITY)
    elif fmt == "avif":
        img.save(buf, format="AVIF", quality=AVIF_QUALITY)
    else:
        raise ValueError(f"format de tuile inconnu : {fmt}")
    return buf.getvalue()


def _write_atomic(path, data):
    path.parent.mkdir(parents=True, exist_ok=True)
    tmp = path.with_name(path.name + f".{os.getpid()}.part")
    try:
        tmp.write_bytes(data)
        os.replace(tmp, path)
    except OSError as e:
        logger.debug(f"Écriture de tuile impossible ({path}) : {e}")
        tmp.unlink(missing_ok=True)


def zoom_supported(z, scale=1):
    """Vrai si le zoom demandé est dans la plage rendue par le serveur."""
    max_z = TILE_MAX_NATIVE_Z - (1 if scale == 2 else 0)
    return TILE_MIN_Z <= z <= max_z


def get_tile(output_dir, layer, z, x, y, scale=1, fmt="png"):
    """Tuile encodée (bytes) depuis le cache, rendue au besoin.

    Retourne None quand aucune donnée n'intersecte la tuile ; un marqueur
    `.empty` mémorise ce cas pour ne jamais le recalculer. Une tuile en cache
    est réutilisée tant qu'aucune dalle contributrice n'est plus récente
    qu'elle : une dalle régénérée n'invalide que ses propres tuiles.
    """
    cache = tile_cache_path(output_dir, layer, z, x, y, scale, fmt)
    empty = cache.with_suffix(".empty")
    sources = _contributing(output_dir, layer, z, x, y, scale)
    if not sources:
        if not empty.exists():
            _write_atomic(empty, b"")
        return None

    newest = 0.0
    for src in sources:
        m = src.mtime()
        if m and m > newest:
            newest = m
    try:
        if cache.stat().st_mtime >= newest:
            return cache.read_bytes()
    except OSError:
        pass

    img = render_tile(output_dir, layer, z, x, y, scale)
    if img is None:
        _write_atomic(empty, b"")
        return None
    data = _encode(img, fmt)
    _write_atomic(cache, data)
    empty.unlink(missing_ok=True)
    return data


def empty_marker_exists(output_dir, layer, z, x, y, scale=1, fmt="png"):
    """Vrai si la tuile est déjà connue comme vide (négatif mémorisé)."""
    return tile_cache_path(output_dir, layer, z, x, y,
                           scale, fmt).with_suffix(".empty").exists()


def cached_tile(output_dir, layer, z, x, y, scale=1, fmt="png"):
    """Tuile depuis le cache UNIQUEMENT — jamais de rendu (mode cache seule).

    Retourne ``(data, état)`` :
    - ``"fresh"``  : tuile en cache et plus récente que ses dalles sources ;
    - ``"empty"``  : aucune source n'intersecte (marqueur `.empty` posé) ;
    - ``"pending"``: donnée présente mais tuile absente ou périmée — à rendre
      ailleurs (tâche de fond, pré-chauffage), pas au fil de la navigation.
    """
    cache = tile_cache_path(output_dir, layer, z, x, y, scale, fmt)
    empty = cache.with_suffix(".empty")
    sources = _contributing(output_dir, layer, z, x, y, scale)
    if not sources:
        if not empty.exists():
            _write_atomic(empty, b"")
        return None, "empty"
    if empty.exists():
        return None, "empty"
    newest = 0.0
    for src in sources:
        m = src.mtime()
        if m and m > newest:
            newest = m
    try:
        if cache.stat().st_mtime >= newest:
            return cache.read_bytes(), "fresh"
    except OSError:
        pass
    return None, "pending"


def transparent_tile(scale=1, fmt="png"):
    """Tuile entièrement transparente (zone sans donnée)."""
    from PIL import Image
    size = TILE_SIZE * scale
    return _encode(Image.new("RGBA", (size, size), (0, 0, 0, 0)), fmt)


# ---------------------------------------------------------------------------
# Pré-chauffage
# ---------------------------------------------------------------------------

def tiles_in_bounds(bounds_wgs84, z):
    """Indices (x, y) des tuiles du niveau z couvrant une bbox WGS84."""
    west, south, east, north = bounds_wgs84
    n = 2 ** z

    def _x(lon):
        return int((lon + 180.0) / 360.0 * n)

    def _y(lat):
        lat = max(-85.0511, min(85.0511, lat))
        rad = math.radians(lat)
        return int((1.0 - math.log(math.tan(rad) + 1 / math.cos(rad)) / math.pi) / 2.0 * n)

    x0, x1 = sorted((_x(west), _x(east)))
    y0, y1 = sorted((_y(north), _y(south)))
    return [(x, y)
            for x in range(max(0, x0), min(n - 1, x1) + 1)
            for y in range(max(0, y0), min(n - 1, y1) + 1)]


def warm(output_dir, layers, z_min, z_max, bounds_wgs84=None, scale=1,
         fmt="png", limit=20000):
    """Pré-calcule les tuiles d'une emprise (retourne le compte rendu)."""
    if bounds_wgs84 is None:
        bounds_wgs84 = grid_bounds_wgs84(output_dir)
    if bounds_wgs84 is None:
        return {"rendues": 0, "vides": 0, "limite": False}
    done = empty = 0
    for layer in layers:
        for z in range(max(TILE_MIN_Z, z_min), min(TILE_MAX_NATIVE_Z, z_max) + 1):
            if not zoom_supported(z, scale):
                continue
            for x, y in tiles_in_bounds(bounds_wgs84, z):
                if done + empty >= limit:
                    return {"rendues": done, "vides": empty, "limite": True}
                if get_tile(output_dir, layer, z, x, y, scale, fmt) is None:
                    empty += 1
                else:
                    done += 1
    return {"rendues": done, "vides": empty, "limite": False}