Le code

def analyser(jeu: JeuCouverture, zone=None, params: Parametres | None = None, progression=None) -> Analyse:
    params = params or Parametres()
    durees = {}

    def etape(nom, debut, fin):
        """Ramène la progression d'une étape dans son intervalle global."""
        if progression is None:
            return None
        return lambda msg, f: progression(msg, debut + (fin - debut) * f)

    t = time.time()
    idx = jeu.selection(zone)
    emprise = shapely.bounds(zone) if zone is not None else jeu.emprise
    res = choisir_resolution(emprise, jeu.resolution_native)
    grille = Grille.depuis_emprise(emprise, res)
    L = rasteriser_niveaux(jeu.parts[idx], jeu.niveau[idx], grille, zone, etape("raster", 0.0, 0.2))
    if params.ignorer_exterieur:
        if progression:
            progression("Repérage de l'extérieur (mer, bords)", 0.2)
        L = marquer_exterieur(L)
    durees["rasterisation"] = round(time.time() - t, 2)

    t = time.time()
    utile = float((L <= 3).sum()) * res * res
    pas = params.pas or taille_automatique(utile, res)
    maillage = Maillage.couvrant(grille.emprise, pas)
    if progression:
        progression(f"Maillage hexagonal ({pas:.0f} m)", 0.25)
    durees["maillage"] = round(time.time() - t, 2)

    t = time.time()
    R = calculer_descripteurs(L, grille, maillage, jeu.binaire, etape("desc", 0.25, 0.75))
    durees["descripteurs"] = round(time.time() - t, 2)

    t = time.time()
    if progression:
        progression("Score (ECOD, Isolation Forest, kNN)", 0.78)
    ignorer = ("transitions_brutales",) if jeu.binaire else ()
    S = scorer(R.X, R.a_scorer, ignorer=ignorer)
    durees["score"] = round(time.time() - t, 2)

    t = time.time()
    if progression:
        progression("Extraction des objets d'anomalie", 0.8)
    O = Objets(L, R, grille, maillage, jeu.binaire)
    durees["objets"] = round(time.time() - t, 2)

    a = Analyse(jeu, zone, params, grille, L, maillage, R, S, O, durees)

    if params.controle:
        t = time.time()
        if progression:
            progression("Contrôle : injection d'anomalies synthétiques", 0.82)
        L2, poses = ctrl.injecter(L, grille, maillage, R, S)
        R2 = calculer_descripteurs(L2, grille, maillage, jeu.binaire, etape("ctrl", 0.82, 0.97))
        S2 = scorer(R2.X, R2.a_scorer, ignorer=ignorer)
        O2 = Objets(L2, R2, grille, maillage, jeu.binaire)
        a.controle = ctrl.evaluer(poses, R2, S2, O2, params.seuil)
        del L2, R2, O2
        durees["controle"] = round(time.time() - t, 2)

    if progression:
        progression("Terminé", 1.0)
    return a

def lire_couverture(chemin: str | Path, progression=None) -> JeuCouverture:
    """Lit un GeoPackage Arcep et l'éclate en polygones simples."""
    chemin = Path(chemin)
    if progression:
        progression("Lecture du fichier", 0.05)

    df = pyogrio.read_dataframe(chemin)
    if df.crs is None:
        raise ValueError("Le fichier n'a pas de système de coordonnées déclaré.")
    crs = CRS.from_user_input(df.crs)

    # On travaille toujours dans un CRS métrique. Les fichiers Arcep le sont
    # déjà (Lambert-93 pour la métropole, UTM pour l'outre-mer), mais on se
    # protège au cas où quelqu'un donnerait un fichier en degrés.
    if crs.is_geographic:
        crs_metrique = df.estimate_utm_crs()
        df = df.to_crs(crs_metrique)
        crs = CRS.from_user_input(crs_metrique)

    if progression:
        progression("Éclatement des multipolygones", 0.5)

    if "niveau" in df.columns:
        codes = df["niveau"].fillna("").astype(str).str.upper().map(NIVEAUX)
    else:
        codes = np.full(len(df), np.nan)
    binaire = bool(np.all(np.isnan(codes.to_numpy(dtype=float))))
    if binaire:
        # Fichier sans niveaux (5G par exemple) : tout ce qui est couvert est
        # rangé dans le niveau le plus haut, faute de mieux.
        codes = np.full(len(df), 3)
    codes = np.nan_to_num(np.asarray(codes, dtype=float), nan=1).astype(np.uint8)

    parts, idx = shapely.get_parts(np.asarray(df.geometry.values), return_index=True)
    # On ne garde que les polygones (un fichier mal formé peut contenir des
    # lignes ou des points perdus, qui n'ont pas de sens pour une couverture).
    garder = shapely.get_type_id(parts) == 3
    parts, idx = parts[garder], idx[garder]
    niveau = codes[idx]

    attributs = {}
    for col in ("operateur_commercial", "techno", "usage", "date", "dept"):
        if col in df.columns:
            attributs[col] = str(df[col].iloc[0])
    del df  # on ne garde que les morceaux, la table d'origine prend de la place

    if progression:
        progression("Index spatial", 0.85)
    arbre = shapely.STRtree(parts)

    return JeuCouverture(
        chemin=chemin,
        crs=crs,
        parts=parts,
        niveau=niveau,
        binaire=binaire,
        resolution_native=_estimer_resolution(parts),
        attributs=attributs,
        arbre=arbre,
    )

def _estimer_resolution(parts: np.ndarray, n_max: int = 3000) -> float:
    """Estime le pas du raster sous-jacent à partir des segments des contours.

    On prend le premier quartile des longueurs de segments horizontaux ou
    verticaux. Sur les fichiers de La Réunion on retrouve exactement 10 m. Si
    le fichier n'est pas en escalier (vrais polygones vectoriels), on obtient
    une valeur plus grande et on la borne à 25 m pour rester raisonnable.
    """
    echantillon = parts[: min(n_max, len(parts))]
    anneaux = shapely.get_exterior_ring(echantillon)
    coords, idx = shapely.get_coordinates(anneaux, return_index=True)
    d = np.diff(coords, axis=0)[idx[1:] == idx[:-1]]
    longueur = np.hypot(d[:, 0], d[:, 1])
    aligne = (np.abs(d[:, 0]) < 1e-6) | (np.abs(d[:, 1]) < 1e-6)
    longueur = longueur[aligne & (longueur > 0.5)]
    if len(longueur) < 50:
        return 25.0
    pas = float(np.quantile(longueur, 0.25))
    return float(np.clip(round(pas, 1), 1.0, 25.0))

def rasteriser_parite(geoms: np.ndarray, grille: Grille, bloc: int = 1024) -> np.ndarray:
    """Masque booléen H x W des pixels dont le centre est dans les polygones."""
    masque = np.zeros((grille.H, grille.W), dtype=bool)
    if len(geoms) == 0:
        return masque

    coords, idx = _anneaux(geoms)
    # Passage en coordonnées pixel : colonnes vers la droite, lignes vers le bas.
    px = (coords[:, 0] - grille.x0) / grille.res
    py = (grille.y1 - coords[:, 1]) / grille.res

    meme = idx[1:] == idx[:-1]  # segments internes à un même anneau
    xa, ya = px[:-1][meme], py[:-1][meme]
    xb, yb = px[1:][meme], py[1:][meme]

    # Un segment coupe la ligne de centres r + 0.5 si lo <= r + 0.5 < hi.
    lo = np.minimum(ya, yb)
    hi = np.maximum(ya, yb)
    r0 = np.ceil(lo - 0.5).astype(np.int64)
    r1 = np.ceil(hi - 0.5).astype(np.int64)
    utile = r1 > r0  # les segments horizontaux ne coupent aucune ligne
    xa, ya, xb, yb, r0, r1 = xa[utile], ya[utile], xb[utile], yb[utile], r0[utile], r1[utile]
    pente = (xb - xa) / (yb - ya)

    W1 = grille.W + 1
    for b0 in range(0, grille.H, bloc):
        b1 = min(b0 + bloc, grille.H)
        sel = (r0 < b1) & (r1 > b0)
        if not sel.any():
            continue
        s0 = np.maximum(r0[sel], b0)
        s1 = np.minimum(r1[sel], b1)
        n = s1 - s0
        # Développement : une entrée par couple (segment, ligne traversée).
        rep = np.repeat(np.arange(len(n)), n)
        decal = np.arange(n.sum()) - np.repeat(np.cumsum(n) - n, n)
        r = s0[rep] + decal
        x = xa[sel][rep] + (r + 0.5 - ya[sel][rep]) * pente[sel][rep]
        # Le pixel c est à droite du croisement si c + 0.5 > x.
        c = np.clip(np.floor(x - 0.5).astype(np.int64) + 1, 0, grille.W)
        compte = np.bincount((r - b0) * W1 + c, minlength=(b1 - b0) * W1)
        bascules = compte.astype(np.uint8).reshape(b1 - b0, W1)
        # La parité d'une somme cumulée en uint8 reste juste malgré le
        # débordement (256 est pair), ce qui évite des tableaux en int64.
        masque[b0:b1] = (np.cumsum(bascules, axis=1, dtype=np.uint8)[:, : grille.W] & 1).astype(bool)
    return masque

def marquer_exterieur(L: np.ndarray) -> np.ndarray:
    """Passe en EXTERIEUR les zones sans couverture reliées au bord.

    Première version du scoring : toute la côte de La Réunion ressortait en
    "transition brutale", parce que la très bonne couverture s'arrête net sur
    le trait de côte. Ce n'est pas une anomalie, c'est le découpage à la mer
    fait lors de la publication. On considère donc qu'une zone sans couverture
    qui touche le bord de la sélection (ou la sélection elle-même) est de
    l'extérieur, pas un trou dans la couverture.

    Limite qu'on assume : une grande zone non couverte à l'intérieur des
    terres mais qui touche le bord du rectangle sélectionné sera aussi
    considérée comme extérieure. Il suffit de sélectionner plus large.
    """
    vide = L == 0
    lab, n = ndimage.label(vide)
    if n == 0:
        return L
    bord = np.zeros(n + 1, dtype=bool)
    for tranche in (lab[0], lab[-1], lab[:, 0], lab[:, -1]):
        bord[np.unique(tranche)] = True
    # Contact avec le hors-zone (cas d'une sélection au lasso).
    hz = L == HORS_ZONE
    if hz.any():
        contact = ndimage.binary_dilation(hz) & vide
        bord[np.unique(lab[contact])] = True
    bord[0] = False
    L = L.copy()
    L[bord[lab]] = EXTERIEUR
    return L

@dataclass
class Maillage:
    """Grille hexagonale à sommet plat, indexée de façon dense."""

    pas: float  # distance entre centres voisins (m)
    ox: float  # origine locale (m), pour garder des petits nombres
    oy: float
    qmin: int
    rmin: int
    nq: int
    nr: int

    @property
    def rayon(self) -> float:
        """Rayon du cercle circonscrit (centre vers sommet)."""
        return self.pas / RAC3

    @property
    def aire(self) -> float:
        return self.pas**2 * RAC3 / 2

    @property
    def n(self) -> int:
        return self.nq * self.nr

    @classmethod
    def couvrant(cls, emprise, pas: float) -> Maillage:
        minx, miny, maxx, maxy = emprise
        m = cls(pas, minx, miny, 0, 0, 1, 1)
        coins_x = np.array([minx, maxx, minx, maxx])
        coins_y = np.array([miny, miny, maxy, maxy])
        q, r = m._axial_brut(coins_x, coins_y)
        qmin, qmax = int(np.floor(q.min())) - 2, int(np.ceil(q.max())) + 2
        rmin, rmax = int(np.floor(r.min())) - 2, int(np.ceil(r.max())) + 2
        return cls(pas, minx, miny, qmin, rmin, qmax - qmin + 1, rmax - rmin + 1)

    def _axial_brut(self, x, y):
        a = self.rayon
        x = np.asarray(x, dtype=np.float64) - self.ox
        y = np.asarray(y, dtype=np.float64) - self.oy
        q = (2.0 / 3.0) * x / a
        r = (-x / 3.0 + RAC3 / 3.0 * y) / a
        return q, r

    def indice(self, x, y) -> np.ndarray:
        """Indice dense de la maille contenant chaque point (arrondi cubique)."""
        q, r = self._axial_brut(x, y)
        s = -q - r
        rq, rr, rs = np.round(q), np.round(r), np.round(s)
        dq, dr, ds = np.abs(rq - q), np.abs(rr - r), np.abs(rs - s)
        corr_q = (dq > dr) & (dq > ds)
        corr_r = ~corr_q & (dr > ds)
        rq = np.where(corr_q, -rr - rs, rq)
        rr = np.where(corr_r, -rq - rs, rr)
        iq = rq.astype(np.int64) - self.qmin
        ir = rr.astype(np.int64) - self.rmin
        return iq * self.nr + ir

    def qr(self, ids):
        ids = np.asarray(ids)
        return ids // self.nr + self.qmin, ids % self.nr + self.rmin

    def centres(self, ids):
        q, r = self.qr(ids)
        a = self.rayon
        x = a * 1.5 * q + self.ox
        y = a * RAC3 * (r + q / 2.0) + self.oy
        return x, y

    def sommets(self, ids) -> np.ndarray:
        """Tableau (n, 7, 2) des sommets de chaque hexagone, anneau fermé."""
        cx, cy = self.centres(ids)
        ang = np.deg2rad(np.arange(0, 420, 60))
        a = self.rayon
        return np.stack([cx[:, None] + a * np.cos(ang), cy[:, None] + a * np.sin(ang)], axis=-1)

    def voisins(self, ids) -> np.ndarray:
        """Tableau (n, 6) des indices voisins, -1 si hors de la grille."""
        q, r = self.qr(ids)
        qn = q[:, None] + DIRECTIONS[:, 0]
        rn = r[:, None] + DIRECTIONS[:, 1]
        iq, ir = qn - self.qmin, rn - self.rmin
        ok = (iq >= 0) & (iq < self.nq) & (ir >= 0) & (ir < self.nr)
        return np.where(ok, iq * self.nr + ir, -1)

def taille_automatique(surface_utile_m2: float, res: float, cible: int = 12000, pixels_min: int = 900) -> float:
    """Pas (distance entre centres voisins, en m) arrondi à une valeur ronde."""
    aire_cible = surface_utile_m2 / cible
    aire_min = pixels_min * res * res
    aire = max(aire_cible, aire_min)
    # Pour un hexagone de pas p (= distance entre deux centres voisins),
    # l'aire vaut p^2 * sqrt(3) / 2.
    pas = math.sqrt(2 * aire / RAC3)
    for t in TAILLES_RONDES:
        if t >= pas:
            return float(t)
    return float(TAILLES_RONDES[-1])

def calculer_descripteurs(L: np.ndarray, grille: Grille, maillage: Maillage, binaire: bool = False, progression=None, bloc: int = 256) -> ResultatDescripteurs:
    H, W, res = grille.H, grille.W, grille.res
    n = maillage.n
    drapeaux = np.zeros((H, W), dtype=np.uint8)
    Lp = np.pad(L, 1, constant_values=HORS_ZONE)

    comptes = np.zeros(n * 4, dtype=np.int64)
    paires = np.zeros(n, dtype=np.int64)
    diff = np.zeros(n, dtype=np.int64)
    brut = np.zeros(n, dtype=np.int64)
    rect_long = np.zeros(n, dtype=np.float64)
    rect_max = np.zeros(n, dtype=np.float64)
    seuil_run = max(4, math.ceil(LONGUEUR_RECTILIGNE_M / res))
    hist_runs = np.zeros(1, dtype=np.int64)
    runs = {k: [] for k in ("sens", "r", "c", "lg", "niv_a", "niv_b")}

    def ajouter_runs(rows, cols, longueurs, sens):
        """Accumule les longues frontières droites dans les mailles."""
        nonlocal hist_runs
        h = np.bincount(longueurs)
        if len(h) > len(hist_runs):
            h[: len(hist_runs)] += hist_runs
            hist_runs = h
        else:
            hist_runs[: len(h)] += h
        garder = longueurs >= seuil_run
        rows, cols, longueurs = rows[garder], cols[garder], longueurs[garder]
        if not len(longueurs):
            return
        if sens == "h":  # frontière horizontale : la suite avance en colonnes
            x = grille.x0 + (cols + longueurs / 2) * res
            y = grille.y1 - (rows + 1) * res
            pr = np.repeat(rows, longueurs)
            pc = np.repeat(cols, longueurs) + np.arange(longueurs.sum()) - np.repeat(np.cumsum(longueurs) - longueurs, longueurs)
        else:  # frontière verticale : la suite avance en lignes
            x = grille.x0 + (cols + 1) * res
            y = grille.y1 - (rows + longueurs / 2) * res
            pc = np.repeat(cols, longueurs)
            pr = np.repeat(rows, longueurs) + np.arange(longueurs.sum()) - np.repeat(np.cumsum(longueurs) - longueurs, longueurs)
        # Une frontière droite qui longe une bande d'un pixel est le bord de
        # cette bande : on l'attribue à la bande fine, pas aux bords
        # rectilignes, sinon le même défaut serait compté deux fois.
        cote_b = (drapeaux[np.minimum(pr + 1, H - 1), pc] if sens == "h" else drapeaux[pr, np.minimum(pc + 1, W - 1)]) & BIT_BANDE
        fin = ((drapeaux[pr, pc] & BIT_BANDE) | cote_b) > 0
        debut = np.concatenate([[0], np.cumsum(longueurs)[:-1]])
        part_fine = np.add.reduceat(fin.astype(np.int64), debut) / longueurs
        garder = part_fine <= 0.5
        if not garder.all():
            garde_px = np.repeat(garder, longueurs)
            rows, cols, longueurs, x, y = rows[garder], cols[garder], longueurs[garder], x[garder], y[garder]
            pr, pc = pr[garde_px], pc[garde_px]
            if not len(longueurs):
                return
        ids = maillage.indice(x, y)
        # Niveaux de part et d'autre, lus au milieu de la frontière :
        # nord / sud pour une frontière horizontale, ouest / est sinon.
        if sens == "h":
            mr, mc = rows, cols + longueurs // 2
            na, nb = L[mr, mc], L[np.minimum(mr + 1, H - 1), mc]
        else:
            mr, mc = rows + longueurs // 2, cols
            na, nb = L[mr, mc], L[mr, np.minimum(mc + 1, W - 1)]
        for k, v in zip(("sens", "r", "c", "lg", "niv_a", "niv_b"), (np.full(len(rows), sens), rows, cols, longueurs, na, nb), strict=True):
            runs[k].append(np.asarray(v))
        np.add.at(rect_long, ids, longueurs * res)
        np.maximum.at(rect_max, ids, longueurs * res)
        drapeaux[pr, pc] |= BIT_RECTILIGNE

    # Passe 1 : composantes connexes ----------------------------------------
    # Faite en premier parce que les passes suivantes ont besoin de savoir
    # quels pixels appartiennent à des micro-zones : un pixel isolé est
    # forcément "fin" et forcément en transition avec ses voisins, on ne veut
    # pas le compter trois fois sous trois noms différents.
    micro_px = max(1, round(SURFACE_MICRO_M2 / res**2))
    enclave_px = max(micro_px, round(SURFACE_ENCLAVE_MAX_M2 / res**2))
    n_micro = np.zeros(n, dtype=np.float64)
    isolees = np.zeros(n, dtype=np.float64)
    n_trous = np.zeros(n, dtype=np.int64)
    n_ilots = np.zeros(n, dtype=np.int64)
    enclaves = {k: [] for k in ("niveau", "entoure", "taille", "graine", "r0", "r1", "c0", "c1")}
    n_zones = 0
    niveaux_presents = (0, 3) if binaire else (0, 1, 2, 3)
    for i, k in enumerate(niveaux_presents):
        if progression:
            progression(f"Zones isolées (niveau {k})", 0.4 * i / len(niveaux_presents))
        lab, nlab = ndimage.label(L == k)
        if nlab == 0:
            continue
        n_zones += nlab
        taille = np.bincount(lab.ravel(), minlength=nlab + 1)
        petit = taille <= enclave_px
        petit[0] = False

        # Classes voisines de chaque petite composante, codées en bits
        # (bit 4 = bord de la sélection ou extérieur, ce qui interdit d'être
        # une enclave : une zone qui touche la mer n'est pas "entourée").
        # (On travaille par tranches décalées plutôt que sur une copie
        # agrandie des étiquettes : ça économise 150 Mo sur La Réunion, ce qui
        # compte sur un petit serveur.)
        bits = np.zeros(nlab + 1, dtype=np.uint8)
        for A, V in (
            (lab[:, :-1], L[:, 1:]),
            (lab[:, 1:], L[:, :-1]),
            (lab[:-1, :], L[1:, :]),
            (lab[1:, :], L[:-1, :]),
        ):
            m = (A > 0) & (V != k)
            m &= petit[A]
            vois = V[m]
            np.bitwise_or.at(bits, A[m], np.where(vois > 3, 16, 1 << np.minimum(vois, 3)).astype(np.uint8))
            del m, vois
        for bord in (lab[0], lab[-1], lab[:, 0], lab[:, -1]):
            bits[np.unique(bord)] |= 16

        unique = (bits != 0) & ((bits & (bits - 1)) == 0) & (bits < 16)
        entoure = np.where(unique, np.log2(np.maximum(bits, 1)).astype(np.int16), -1)
        contraste = np.where(unique, entoure - k, 0)
        # On ne garde que les enclaves avec un saut d'au moins deux niveaux.
        # Une poche de "bonne" dans de la "très bonne" couverture, c'est le
        # dégradé normal du signal, on en trouve partout en montagne.
        # Les micro-zones (4 pixels ou moins) sont comptées à part : une zone
        # isolée, c'est une vraie zone, pas un pixel de bruit.
        micro = petit & (taille <= micro_px)
        enclave = petit & ~micro & unique & (np.abs(contraste) >= 2)

        # Une zone plus étroite que deux pixels en moyenne est une bande, pas
        # une zone isolée : largeur moyenne = surface / plus grande dimension.
        cand = np.flatnonzero(enclave)
        if len(cand):
            pix = np.flatnonzero(enclave[lab.ravel()])
            lb = lab.ravel()[pix]
            pr_, pc_ = pix // W, pix % W
            h0 = np.full(nlab + 1, H)
            h1 = np.zeros(nlab + 1, dtype=np.int64)
            w0 = np.full(nlab + 1, W)
            w1 = np.zeros(nlab + 1, dtype=np.int64)
            np.minimum.at(h0, lb, pr_)
            np.maximum.at(h1, lb, pr_)
            np.minimum.at(w0, lb, pc_)
            np.maximum.at(w1, lb, pc_)
            etendue = np.maximum(h1[cand] - h0[cand], w1[cand] - w0[cand]) + 1
            enclave[cand[taille[cand] / etendue < 2]] = False

        interet = enclave | micro
        idx = np.flatnonzero(interet[lab.ravel()])
        if len(idx) == 0:
            continue
        labs = lab.ravel()[idx]
        r, c = idx // W, idx % W
        cr = np.bincount(labs, weights=r, minlength=nlab + 1)
        cc = np.bincount(labs, weights=c, minlength=nlab + 1)
        comp = np.flatnonzero(interet)
        x = grille.x0 + (cc[comp] / taille[comp] + 0.5) * res
        y = grille.y1 - (cr[comp] / taille[comp] + 0.5) * res
        cell = maillage.indice(x, y)
        np.add.at(n_micro, cell[micro[comp]], 1)
        e = enclave[comp]
        np.add.at(isolees, cell[e], np.abs(contraste[comp][e]) * taille[comp][e])
        np.add.at(n_trous, cell[e & (contraste[comp] > 0)], 1)
        np.add.at(n_ilots, cell[e & (contraste[comp] < 0)], 1)

        # On garde chaque enclave comme un objet : un pixel graine (pour la
        # retrouver plus tard) et sa boîte englobante.
        graine = np.full(nlab + 1, np.iinfo(np.int64).max)
        np.minimum.at(graine, labs, idx)
        r0 = np.full(nlab + 1, H)
        r1 = np.zeros(nlab + 1, dtype=np.int64)
        c0 = np.full(nlab + 1, W)
        c1 = np.zeros(nlab + 1, dtype=np.int64)
        np.minimum.at(r0, labs, r)
        np.maximum.at(r1, labs, r)
        np.minimum.at(c0, labs, c)
        np.maximum.at(c1, labs, c)
        ce = np.flatnonzero(enclave)
        for k_, v in zip(("niveau", "entoure", "taille", "graine", "r0", "r1", "c0", "c1"), (np.full(len(ce), k), entoure[ce], taille[ce], graine[ce], r0[ce], r1[ce], c0[ce], c1[ce]), strict=True):
            enclaves[k_].append(np.asarray(v))

        flat = drapeaux.ravel()
        flat[idx[enclave[labs]]] |= BIT_ENCLAVE
        flat[idx[micro[labs]]] |= BIT_MICRO
        del lab

    micro_p = np.pad((drapeaux & BIT_MICRO) > 0, 1)

    # Passe 2a : bandes fines, sur tout le raster, avant les transitions.
    # Un pixel de bande est différent de ses deux voisins opposés, eux-mêmes
    # identiques : c'est la signature d'un polygone d'un pixel de large.
    for b0 in range(0, H, bloc):
        b1 = min(b0 + bloc, H)
        C = Lp[b0 + 1 : b1 + 1, 1:-1]
        haut, bas = Lp[b0:b1, 1:-1], Lp[b0 + 2 : b1 + 2, 1:-1]
        gauche, droite = Lp[b0 + 1 : b1 + 1, :-2], Lp[b0 + 1 : b1 + 1, 2:]
        ok = (C <= 3) & ~micro_p[b0 + 1 : b1 + 1, 1:-1]
        fine_h = ok & (gauche == droite) & (gauche != C) & (gauche <= 3)
        fine_v = ok & (haut == bas) & (haut != C) & (haut <= 3)
        drapeaux[b0:b1][fine_h | fine_v] |= BIT_BANDE
    # Pixels exclus des transitions : micro-zones, bandes et zones isolées.
    # Chaque défaut n'est compté que sous un seul nom : le bord d'un trou
    # n'est pas une « transition brutale », c'est le trou qui est anormal.
    exclu_p = micro_p | np.pad((drapeaux & (BIT_BANDE | BIT_ENCLAVE)) > 0, 1)

    # Bandes fines significatives : on regroupe les pixels de bande en objets
    # (8-connexité) et on ne garde que ceux d'au moins 80 m. Compter tous les
    # pixels fins noyait une vraie bande dans le bruit de la frontière.
    bande_min = max(8, math.ceil(80 / res))
    lab_b, nb_b = ndimage.label((drapeaux & BIT_BANDE) > 0, structure=np.ones((3, 3), dtype=bool))
    long_bandes = np.zeros(n, dtype=np.float64)
    px_bandes = np.zeros(n, dtype=np.int64)
    if nb_b:
        t_b = np.bincount(lab_b.ravel(), minlength=nb_b + 1)
        pix = np.flatnonzero(lab_b.ravel())
        lb = lab_b.ravel()[pix]
        cr_b = np.bincount(lb, weights=pix // W, minlength=nb_b + 1)
        cc_b = np.bincount(lb, weights=pix % W, minlength=nb_b + 1)
        sig = np.flatnonzero(t_b >= bande_min)
        sig = sig[sig > 0]
        cell_b = maillage.indice(grille.x0 + (cc_b[sig] / t_b[sig] + 0.5) * res, grille.y1 - (cr_b[sig] / t_b[sig] + 0.5) * res)
        np.add.at(long_bandes, cell_b, t_b[sig] * res)
        np.add.at(px_bandes, cell_b, t_b[sig])
    del lab_b

    # Passe 2b : par blocs de lignes ---------------------------------------
    for b0 in range(0, H, bloc):
        b1 = min(b0 + bloc, H)
        if progression:
            progression("Transitions et bandes fines", 0.4 + 0.4 * b0 / H)
        X, Y = _coord_pixels(grille, b0, b1)
        ids = maillage.indice(X, Y)
        C = Lp[b0 + 1 : b1 + 1, 1:-1]
        haut = Lp[b0 : b1, 1:-1]
        bas = Lp[b0 + 2 : b1 + 2, 1:-1]
        gauche = Lp[b0 + 1 : b1 + 1, :-2]
        droite = Lp[b0 + 1 : b1 + 1, 2:]
        ok = C <= 3
        # Pixels de micro-zones ou de bandes et leurs voisins sont exclus des
        # transitions (déjà comptés sous leur propre nom).
        Mp = exclu_p[b0 : b1 + 2]
        net = ok & ~Mp[1:-1, 1:-1]
        nets = (Mp[1:-1, 2:], Mp[2:, 1:-1], Mp[1:-1, :-2], Mp[:-2, 1:-1])

        comptes += np.bincount((ids[ok] * 4 + C[ok]).ravel(), minlength=n * 4)

        Ci = C.astype(np.int16)
        brut_px = np.zeros(C.shape, dtype=bool)
        for voisin, mv in ((droite, nets[0]), (bas, nets[1])):  # une fois par paire
            v = ok & (voisin <= 3)
            d = np.abs(Ci - voisin.astype(np.int16))
            paires += np.bincount(ids[v], minlength=n)
            diff += np.bincount(ids[v & (d > 0)], minlength=n)
            b = v & (d >= 2) & net & ~mv
            brut += np.bincount(ids[b], minlength=n)
            brut_px |= b
        for voisin, mv in ((gauche, nets[2]), (haut, nets[3])):
            v = ok & (voisin <= 3)
            brut_px |= v & (np.abs(Ci - voisin.astype(np.int16)) >= 2) & net & ~mv

        if not binaire:
            drapeaux[b0:b1][brut_px] |= BIT_BRUTAL


    # Passe 3 : frontières droites. Faite après la passe 2 pour que toutes
    # les bandes fines soient déjà repérées. Horizontales par blocs de lignes,
    # verticales par blocs de colonnes.
    for b0 in range(0, H, bloc):
        b1 = min(b0 + bloc, H)
        if progression:
            progression("Frontières rectilignes", 0.8 + 0.05 * b0 / H)
        C = L[b0:b1]
        bas = Lp[b0 + 2 : b1 + 2, 1:-1]
        E = (C <= 3) & (bas <= 3) & (bas != C)
        rs, cs, lg = _runs(E)
        ajouter_runs(rs + b0, cs, lg, "h")

    for c0 in range(0, W - 1, bloc):
        c1 = min(c0 + bloc, W - 1)
        if progression:
            progression("Frontières rectilignes", 0.85 + 0.05 * c0 / W)
        A = L[:, c0:c1]
        B = L[:, c0 + 1 : c1 + 1]
        E = (A <= 3) & (B <= 3) & (A != B)
        rs, cs, lg = _runs(np.ascontiguousarray(E.T))  # lignes de E.T = colonnes
        ajouter_runs(cs, rs + c0, lg, "v")

    # Assemblage par maille -------------------------------------------------
    if progression:
        progression("Assemblage par maille", 0.92)
    comptes = comptes.reshape(n, 4)
    n_px = comptes.sum(axis=1)
    attendu = maillage.aire / res**2
    actives = np.flatnonzero(n_px >= 0.3 * attendu)

    npx = n_px[actives].astype(np.float64)
    fr = comptes[actives] / npx[:, None]
    qualite = (fr * np.arange(4)).sum(axis=1) / 3.0
    couverte = fr[:, 0] < 1.0

    # Voisinage : écart entre la qualité de la maille et celle de ses voisines.
    pos = np.full(n, -1)
    pos[actives] = np.arange(len(actives))
    vois = maillage.voisins(actives)
    vpos = np.where(vois >= 0, pos[np.maximum(vois, 0)], -1)
    valide = vpos >= 0
    qv = np.where(valide, qualite[np.maximum(vpos, 0)], 0.0)
    nb = valide.sum(axis=1)
    moy_vois = np.where(nb > 0, qv.sum(axis=1) / np.maximum(nb, 1), qualite)
    ecart = np.abs(qualite - moy_vois)
    couverte_vois = (np.where(valide, couverte[np.maximum(vpos, 0)], False)).any(axis=1)

    d_ = diff[actives].astype(np.float64)
    X = {
        "transitions_brutales": np.zeros(len(actives)) if binaire else brut[actives] / (d_ + LISSAGE_TRANSITIONS),
        "densite_frontieres": d_ / np.maximum(paires[actives], 1),
        # Micro-zones rapportées à la longueur de frontière (en km) : une zone
        # de relief a naturellement beaucoup de frontières, ce qu'on cherche
        # c'est une frontière anormalement "bruitée" en confettis.
        "micro_zones": n_micro[actives] / (d_ * res / 1000.0 + 0.5),
        "zones_isolees": isolees[actives] / npx,
        "bandes_fines": long_bandes[actives] / maillage.pas,
        "bords_rectilignes": rect_long[actives] / maillage.pas,
        "ecart_voisinage": ecart,
    }
    detail = {
        "paires": paires[actives],
        "frontieres": diff[actives],
        "brutales": brut[actives],
        "n_micro": n_micro[actives].astype(np.int64),
        "n_trous": n_trous[actives],
        "n_ilots": n_ilots[actives],
        "px_fins": px_bandes[actives],
        "rect_long_m": rect_long[actives],
        "rect_max_m": rect_max[actives],
        "qualite_voisins": moy_vois,
        "n_voisins": nb,
    }
    # On ne compare au modèle que les mailles où il se passe quelque chose :
    # couvertes, ou voisines d'une maille couverte (un trou franc en plein
    # milieu d'une zone couverte doit pouvoir remonter). La mer et les zones
    # vides loin de tout restent hors du calcul, sinon elles écraseraient la
    # distribution "normale".
    a_scorer = couverte | couverte_vois

    if progression:
        progression("Descripteurs terminés", 1.0)
    cat = lambda d: {k: (np.concatenate(v) if v else np.array([], dtype=np.int64)) for k, v in d.items()}  # noqa: E731
    return ResultatDescripteurs(actives, n_px[actives], fr, qualite, X, detail, a_scorer, drapeaux, cat(runs), hist_runs, cat(enclaves), n_zones)

def _runs(E: np.ndarray):
    """Débuts et longueurs des suites de True sur chaque ligne de E."""
    Ep = np.zeros((E.shape[0], E.shape[1] + 2), dtype=np.int8)
    Ep[:, 1:-1] = E
    d = np.diff(Ep, axis=1)
    rs, cs = np.nonzero(d == 1)
    re, ce = np.nonzero(d == -1)
    return rs, cs, ce - cs

def _ecod_droite(Xs: np.ndarray):
    """Contributions ECOD sur la queue droite et percentiles empiriques."""
    n, d = Xs.shape
    contrib = np.zeros_like(Xs)
    pct = np.zeros_like(Xs)
    for j in range(d):
        v = Xs[:, j]
        tri = np.sort(v)
        sup_ou_egal = n - np.searchsorted(tri, v, side="left")
        strict_inf = np.searchsorted(tri, v, side="left")
        contrib[:, j] = -np.log(sup_ou_egal / n)
        pct[:, j] = strict_inf / n
    return contrib, pct

def _transformer(Xs: np.ndarray) -> np.ndarray:
    """Mise à l'échelle pour les modèles à distance (IF, kNN).

    Les descripteurs sont très asymétriques (beaucoup de zéros, une longue
    queue). On divise par le 90e centile des valeurs non nulles puis on passe
    en log(1 + x), pour qu'un descripteur ne domine pas les distances juste
    à cause de son unité.
    """
    Z = np.empty_like(Xs)
    for j in range(Xs.shape[1]):
        v = Xs[:, j]
        pos = v[v > 0]
        echelle = np.quantile(pos, 0.9) if len(pos) else 1.0
        Z[:, j] = np.log1p(v / max(echelle, 1e-12))
    return Z

def scorer(X: dict, a_scorer: np.ndarray, graine: int = 42, ignorer: tuple = ()) -> Scores:
    noms = [k for k in DESCRIPTEURS if k not in ignorer]
    Xall = np.column_stack([X[k] for k in noms]).astype(np.float64)
    n = len(Xall)
    sel = np.flatnonzero(a_scorer)
    Xs = Xall[sel]

    proba = np.zeros(n)
    pm = np.zeros((n, 3))
    contrib = np.zeros((n, len(noms)))
    pct = np.zeros((n, len(noms)))
    types = np.full(n, "", dtype=object)
    ordre = np.zeros(n)
    if len(sel) < 30:
        return Scores(proba, pm, np.zeros(n, dtype=int), contrib, pct, types, noms, ordre)

    c, p = _ecod_droite(Xs)
    ctx = [i for i, k in enumerate(noms) if k in CONTEXTE]
    c[:, ctx] = 0.0
    contrib[sel], pct[sel] = c, p
    s_ecod = c.sum(axis=1)

    Z = _transformer(Xs)
    foret = IsolationForest(n_estimators=300, max_samples=min(512, len(sel)), random_state=graine)
    s_if = -foret.fit(Z).score_samples(Z)

    k = min(20, len(sel) - 1)
    dist, _ = NearestNeighbors(n_neighbors=k + 1).fit(Z).kneighbors(Z)
    s_knn = dist[:, 1:].mean(axis=1)

    pm[sel] = np.column_stack([_unifier(s_ecod), _unifier(s_if), _unifier(s_knn)])
    proba[sel] = pm[sel].mean(axis=1)
    rangs = [np.argsort(np.argsort(x)) / (len(x) - 1) for x in (s_ecod, s_if, s_knn)]
    ordre[sel] = np.mean(rangs, axis=0)
    consensus = (pm >= 0.5).sum(axis=1)

    # Type dominant : la famille qui apporte le plus au score ECOD. Si aucun
    # descripteur n'est dans les 5 % les plus extrêmes, c'est la combinaison
    # qui est atypique plutôt qu'un critère précis.
    familles = sorted(set(FAMILLES[k] for k in noms if k not in CONTEXTE))
    par_famille = np.column_stack([contrib[:, [i for i, k in enumerate(noms) if FAMILLES[k] == f and k not in CONTEXTE]].sum(axis=1) for f in familles])
    dominante = np.array(familles, dtype=object)[par_famille.argmax(axis=1)]
    nette = contrib.max(axis=1) >= -np.log(0.05)
    types[sel] = np.where(nette[sel], dominante[sel], TYPE_COMBINE)

    return Scores(proba, pm, consensus, contrib, pct, types, noms, ordre)

def _unifier(s: np.ndarray) -> np.ndarray:
    mu, sigma = s.mean(), s.std()
    if sigma <= 0:
        return np.zeros_like(s)
    return np.clip(erf((s - mu) / (sigma * np.sqrt(2))), 0, 1)

def _composantes(masque: np.ndarray, taille_min: int):
    """Composantes connexes d'un masque : graine, boîte, taille."""
    lab, n = ndimage.label(masque, structure=HUIT)
    if n == 0:
        return None
    taille = np.bincount(lab.ravel(), minlength=n + 1)
    idx = np.flatnonzero(lab.ravel())
    labs = lab.ravel()[idx]
    W = masque.shape[1]
    r, c = idx // W, idx % W
    graine = np.full(n + 1, np.iinfo(np.int64).max)
    np.minimum.at(graine, labs, idx)
    boites = [np.full(n + 1, v, dtype=np.int64) for v in (masque.shape[0], 0, W, 0)]
    np.minimum.at(boites[0], labs, r)
    np.maximum.at(boites[1], labs, r)
    np.minimum.at(boites[2], labs, c)
    np.maximum.at(boites[3], labs, c)
    garder = np.flatnonzero(taille >= taille_min)
    garder = garder[garder > 0]
    return {
        "taille": taille[garder],
        "graine": graine[garder],
        "r0": boites[0][garder],
        "r1": boites[1][garder],
        "c0": boites[2][garder],
        "c1": boites[3][garder],
        "n_total": n,
        "toutes_tailles": taille[1:],
    }

def _rarete_queue(valeurs_toutes: np.ndarray, v: np.ndarray) -> np.ndarray:
    """Part des valeurs de référence supérieures ou égales à v."""
    tri = np.sort(valeurs_toutes)
    return (len(tri) - np.searchsorted(tri, v, side="left")) / max(len(tri), 1)

class Objets:
    """Table des objets d'anomalie d'une analyse."""

    def __init__(self, L, R, grille, maillage, binaire: bool):
        self.L, self.R, self.grille, self.maillage = L, R, grille, maillage
        res = grille.res
        pos = np.full(maillage.n, -1)
        pos[R.ids] = np.arange(len(R.ids))
        self._pos = pos
        lignes = []

        # 1. Segments rectilignes -------------------------------------------
        ru = R.runs
        if len(ru["lg"]):
            total = R.hist_runs.sum()
            cumul = np.cumsum(R.hist_runs[::-1])[::-1]  # nombre de frontières >= longueur
            for i in range(len(ru["lg"])):
                sens, r, c, lg = ru["sens"][i], int(ru["r"][i]), int(ru["c"][i]), int(ru["lg"][i])
                if sens == "h":
                    x, y = grille.x0 + (c + lg / 2) * res, grille.y1 - (r + 1) * res
                    cotes = ("nord", "sud")
                else:
                    x, y = grille.x0 + (c + 1) * res, grille.y1 - (r + lg / 2) * res
                    cotes = ("ouest", "est")
                na, nb = int(ru["niv_a"][i]), int(ru["niv_b"][i])
                rare = cumul[lg] / total
                lignes.append(
                    {
                        "type": "Discontinuité rectiligne",
                        "forme": "segment",
                        "x": x,
                        "y": y,
                        "longueur_m": lg * res,
                        "orientation": "est-ouest" if sens == "h" else "nord-sud",
                        "niv_a": na,
                        "niv_b": nb,
                        "cotes": cotes,
                        "rarete": float(rare),
                        "rare_n": int(cumul[lg]),
                        "rare_total": int(total),
                        "geo": ("segment", sens, r, c, lg),
                        "texte": (
                            f"Frontière parfaitement droite de {lg * res:.0f} m, orientée {('est-ouest' if sens == 'h' else 'nord-sud')}, "
                            f"entre {_niv(na)} au {cotes[0]} et {_niv(nb)} au {cotes[1]}. "
                            f"Sur les {_n(total)} frontières droites de la zone, {_n(cumul[lg])} "
                            f"{_pluriel(int(cumul[lg]), 'est', 'sont')} au moins aussi {_pluriel(int(cumul[lg]), 'longue', 'longues')}."
                        ),
                    }
                )

        # 2. Zones isolées -----------------------------------------------------
        en = R.enclaves
        if len(en["taille"]):
            poids = np.abs(en["entoure"] - en["niveau"]) * en["taille"]
            rare = _rarete_queue(poids, poids) * len(poids) / max(R.n_zones, 1)
            for i in range(len(en["taille"])):
                k, e, t = int(en["niveau"][i]), int(en["entoure"][i]), int(en["taille"][i])
                g = int(en["graine"][i])
                r, c = g // grille.W, g % grille.W
                rc = (en["r0"][i] + en["r1"][i]) / 2 + 0.5
                cc = (en["c0"][i] + en["c1"][i]) / 2 + 0.5
                trou = k < e
                lignes.append(
                    {
                        "type": "Zone isolée",
                        "forme": "trou" if trou else "îlot",
                        "x": grille.x0 + cc * res,
                        "y": grille.y1 - rc * res,
                        "surface_m2": t * res * res,
                        "niv_a": k,
                        "niv_b": e,
                        "rarete": float(rare[i]),
                        "rare_n": int(max(1, round(rare[i] * R.n_zones))),
                        "rare_total": int(R.n_zones),
                        "geo": ("classe", k, r, c, int(en["r0"][i]), int(en["r1"][i]), int(en["c0"][i]), int(en["c1"][i])),
                        "texte": (
                            f"{'Trou' if trou else 'Îlot'} de {_n(t * res * res)} m² en {_niv(k)}, entièrement entouré de {_niv(e)} "
                            f"({abs(e - k)} niveaux d'écart). Sur les {_n(R.n_zones)} zones de la carte, "
                            f"{_n(max(1, round(rare[i] * R.n_zones)))} "
                            f"{_pluriel(max(1, round(rare[i] * R.n_zones)), 'seule est une enclave', 'sont des enclaves')} au moins aussi marquée"
                            f"{_pluriel(max(1, round(rare[i] * R.n_zones)), '', 's')}."
                        ),
                    }
                )

        # 3. Bandes fines et 4. transitions brutales ---------------------------
        min_bande = max(8, int(np.ceil(80 / res)))
        cb = _composantes((R.drapeaux & BIT_BANDE) > 0, min_bande)
        if cb is not None and len(cb["taille"]):
            rare = _rarete_queue(cb["toutes_tailles"], cb["taille"])
            for i in range(len(cb["taille"])):
                lignes.append(self._objet_pixels("Bande fine", "bande", cb, i, rare[i], BIT_BANDE))
        if not binaire:
            cr = _composantes((R.drapeaux & BIT_BRUTAL) > 0, 6)
            if cr is not None and len(cr["taille"]):
                rare = _rarete_queue(cr["toutes_tailles"], cr["taille"])
                for i in range(len(cr["taille"])):
                    lignes.append(self._objet_pixels("Transition brutale", "rupture", cr, i, rare[i], BIT_BRUTAL))

        # Rattachement aux mailles ----------------------------------------------
        self.table = []
        if lignes:
            x = np.array([o["x"] for o in lignes])
            y = np.array([o["y"] for o in lignes])
            cellules = pos[maillage.indice(x, y)]
            for o, cell in zip(lignes, cellules, strict=True):
                if cell >= 0:
                    o["maille"] = int(cell)
                    self.table.append(o)
        for i, o in enumerate(self.table):
            o["id"] = i
        self.par_maille: dict[int, list[int]] = {}
        for o in self.table:
            self.par_maille.setdefault(o["maille"], []).append(o["id"])
        for ids in self.par_maille.values():
            ids.sort(key=lambda j: self.table[j]["rarete"])

    def _objet_pixels(self, type_, forme, comp, i, rare, bit):
        g = self.grille
        res = g.res
        graine = int(comp["graine"][i])
        r, c = graine // g.W, graine % g.W
        r0, r1, c0, c1 = (int(comp[k][i]) for k in ("r0", "r1", "c0", "c1"))
        t = int(comp["taille"][i])
        n_tot = int(comp["n_total"])
        n_sup = max(1, round(rare * n_tot))
        bloc = self.L[r0 : r1 + 1, c0 : c1 + 1]
        dr = self.R.drapeaux[r0 : r1 + 1, c0 : c1 + 1]
        niveaux = np.unique(bloc[((dr & bit) > 0) & (bloc <= 3)])
        etendue = max(r1 - r0 + 1, c1 - c0 + 1) * res
        o = {
            "type": type_,
            "forme": forme,
            "x": g.x0 + ((c0 + c1) / 2 + 0.5) * res,
            "y": g.y1 - ((r0 + r1) / 2 + 0.5) * res,
            "pixels": t,
            "etendue_m": float(etendue),
            "niveaux": [int(v) for v in niveaux],
            "rarete": float(rare),
            "rare_n": int(n_sup),
            "rare_total": int(n_tot),
            "geo": ("drapeau", bit, r, c, r0, r1, c0, c1),
        }
        if type_ == "Bande fine":
            # Niveau de la bande et de ce qui l'entoure, lus à la graine.
            k = int(self.L[r, c])
            gauche = self.L[r, c - 1] if c > 0 else 255
            droite = self.L[r, c + 1] if c + 1 < g.W else 255
            autour = int(gauche if gauche == droite else (self.L[r - 1, c] if r > 0 else 255))
            o["niv_a"], o["niv_b"] = k, autour
            o["texte"] = (
                f"Bande de {t} pixels ({t * res:.0f} m de long environ) et d'un seul pixel ({res:.0f} m) de large, "
                f"en {_niv(k)} au milieu de {_niv(autour)}. Une zone aussi étroite et aussi longue ressemble aux "
                f"« sliver polygons » des superpositions de couches en SIG. Sur les {_n(n_tot)} bandes fines de la zone, "
                f"{_n(n_sup)} {_pluriel(n_sup, 'est', 'sont')} au moins aussi {_pluriel(n_sup, 'longue', 'longues')}."
            )
        else:
            bas, haut = (int(niveaux.min()), int(niveaux.max())) if len(niveaux) else (0, 0)
            o["niv_a"], o["niv_b"] = haut, bas
            o["texte"] = (
                f"Sur environ {t * res / 2:.0f} m de frontière, on passe directement de {_niv(haut)} à {_niv(bas)}, "
                f"sans le palier intermédiaire qu'on attendrait d'un signal qui s'atténue. "
                f"Sur les {_n(n_tot)} ruptures de la zone, {_n(n_sup)} {_pluriel(n_sup, 'est', 'sont')} au moins aussi "
                f"{_pluriel(n_sup, 'étendue', 'étendues')}."
            )
        return o

    # -- géométrie, calculée seulement quand on en a besoin -------------------
    def geometrie(self, oid: int) -> shapely.Geometry:
        g = self.grille
        res = g.res
        geo = self.table[oid]["geo"]
        if geo[0] == "segment":
            _, sens, r, c, lg = geo
            if sens == "h":
                y = g.y1 - (r + 1) * res
                return shapely.LineString([(g.x0 + c * res, y), (g.x0 + (c + lg) * res, y)])
            x = g.x0 + (c + 1) * res
            return shapely.LineString([(x, g.y1 - r * res), (x, g.y1 - (r + lg) * res)])
        quoi, val, r, c, r0, r1, c0, c1 = geo
        bloc = self.L[r0 : r1 + 1, c0 : c1 + 1] == val if quoi == "classe" else (self.R.drapeaux[r0 : r1 + 1, c0 : c1 + 1] & val) > 0
        lab, _ = ndimage.label(bloc, structure=HUIT if quoi == "drapeau" else None)
        m = lab == lab[r - r0, c - c0]
        # Union de rectangles : une boîte par suite horizontale de pixels.
        boites = []
        for i in range(m.shape[0]):
            ligne = np.concatenate([[0], m[i].astype(np.int8), [0]])
            d = np.diff(ligne)
            for a, b in zip(np.flatnonzero(d == 1), np.flatnonzero(d == -1), strict=True):
                y0 = g.y1 - (r0 + i + 1) * res
                boites.append(shapely.box(g.x0 + (c0 + a) * res, y0, g.x0 + (c0 + b) * res, y0 + res))
        return shapely.union_all(boites)

    def resume(self) -> dict:
        types = {}
        for o in self.table:
            types[o["type"]] = types.get(o["type"], 0) + 1
        return types

def injecter(L, grille, maillage, R, S, par_type: int = 10, graine: int = 7):
    """Renvoie une copie de L abîmée et la liste (indice dense, défaut)."""
    rng = np.random.default_rng(graine)
    L2 = L.copy()
    res = grille.res
    poses, interdits = [], set()
    for defaut, masque in candidats(R, S).items():
        ids = rng.permutation(R.ids[masque])
        k = 0
        for cell in ids:
            if k >= par_type:
                break
            # Pas deux défauts dans des mailles voisines : chaque détection
            # doit pouvoir être attribuée à un seul défaut.
            voisins = set(maillage.voisins(np.array([cell]))[0].tolist()) | {int(cell)}
            if voisins & interdits:
                continue
            cx, cy = maillage.centres(np.array([cell]))
            c = int((cx[0] - grille.x0) / res)
            r = int((grille.y1 - cy[0]) / res)
            if defaut == "Trou":
                h, w = max(2, int(80 / res)), max(2, int(120 / res))
                z = L2[r - h // 2 : r + h // 2, c - w // 2 : c + w // 2]
                z[z <= 3] = 0
            elif defaut == "Bande":
                lg = max(4, int(500 / res))
                z = L2[r, c - lg // 2 : c + lg // 2]
                z[z <= 3] = 1
            elif defaut == "Coupure":
                h, w = max(3, int(300 / res)), max(4, int(400 / res))
                z = L2[r - h // 2 : r + h // 2, c - w // 2 : c + w // 2]
                z[z <= 3] = 2
            else:
                d = max(4, int(400 / res)) // 2
                z = L2[r - d : r + d, c - d : c + d]
                z[(z == 1) | (z == 2)] = 0
            poses.append((int(cell), defaut))
            interdits |= voisins
            k += 1
    return L2, poses

def _present(o, L, Lref, R, grille) -> float:
    """Part de l'objet qu'on retrouve à l'identique dans la référence."""
    geo = o["geo"]
    if geo[0] == "segment":
        _, sens, r, c, lg = geo
        if sens == "h":
            a, b = Lref[r, c : c + lg], Lref[min(r + 1, grille.H - 1), c : c + lg]
        else:
            a, b = Lref[r : r + lg, c], Lref[r : r + lg, min(c + 1, grille.W - 1)]
        ok = (a <= 3) & (b <= 3)
        return float(((a != b) & ok).sum() / max(ok.sum(), 1))
    quoi, val, r, c, r0, r1, c0, c1 = geo
    bloc = L[r0 : r1 + 1, c0 : c1 + 1]
    masque = (bloc == val) if quoi == "classe" else ((R.drapeaux[r0 : r1 + 1, c0 : c1 + 1] & val) > 0)
    if not masque.any():
        return 0.0
    ref = Lref[r0 : r1 + 1, c0 : c1 + 1]
    return float((ref[masque] == bloc[masque]).mean())

def enregistrer(a: Analyse, cible: Path):
    # On n'enregistre pas les polygones : ils sont dans le .gpkg et on les
    # relira seulement si quelqu'un demande une sélection plus fine.
    jeu_leger = replace(a.jeu, parts=np.empty(0, dtype=object), niveau=np.empty(0, dtype=np.uint8), arbre=None)
    leger = replace(a, jeu=jeu_leger)
    contenu = {"analyse": leger, "image": a.image(), "version": VERSION}
    cible.parent.mkdir(parents=True, exist_ok=True)
    with gzip.open(cible, "wb", compresslevel=6) as f:
        pickle.dump(contenu, f, protocol=pickle.HIGHEST_PROTOCOL)