gta6/tools/extract/leonida/vectors.py

import networkx as nx
import numpy as np
from scipy import ndimage as ndi
from shapely.geometry import LineString, Polygon
from shapely.ops import unary_union
from skimage.measure import find_contours
from skimage.morphology import closing, disk, skeletonize

from .classify import CLASSES
from .coords import SLICE

ROAD_CLASSES = ("strasse", "weg", "spekulativ")


def _to_world(contour: np.ndarray) -> list[tuple[float, float]]:
    """Konturpunkte (Zeile, Spalte) → Weltkoordinaten in Metern."""
    return [(SLICE["x0"] + c, SLICE["y1"] - r) for r, c in contour]


def mask_polygons(mask: np.ndarray, min_area: float, simplify: float) -> list[Polygon]:
    padded = np.pad(mask, 1)
    outers, holes = [], []
    for contour in find_contours(padded.astype(float), 0.5):
        local = contour - 1
        ring = _to_world(local)
        if len(ring) < 4:
            continue
        poly = Polygon(ring).buffer(0)
        if poly.is_empty or poly.area < min_area or poly.geom_type != "Polygon":
            continue
        (outers if _signed_area(local) > 0 else holes).append(poly)
    result = []
    for outer in outers:
        inner = [h.exterior.coords for h in holes if outer.contains(h)]
        result.append(Polygon(outer.exterior.coords, inner).buffer(0).simplify(simplify))
    return [p for p in result if p.geom_type == "Polygon" and p.area >= min_area]


def _signed_area(contour: np.ndarray) -> float:
    """Shoelace-Summe in Pixelkoordinaten (Zeile, Spalte).

    find_contours liefert Außenkonturen und Löcher mit entgegengesetztem Umlaufsinn;
    das Vorzeichen trennt sie zuverlässig. Eine Stichprobe der Maske am Mittelpunkt der
    Kontur (naheliegende Alternative) versagt bei konzentrischen Formen wie einer Insel
    mit Lagune: der Mittelpunkt der äußeren Kontur liegt dann selbst in der Lagune.
    """
    row, col = contour[:, 0], contour[:, 1]
    return float(0.5 * np.sum(col * np.roll(row, -1) - np.roll(col, -1) * row))


def land_polygons(labels: np.ndarray) -> list[Polygon]:
    land = labels != CLASSES["wasser"]
    land = ndi.binary_closing(ndi.binary_opening(land, disk(2)), disk(2))
    return mask_polygons(land, min_area=50, simplify=1.0)


def road_mask(labels: np.ndarray) -> np.ndarray:
    mask = np.zeros(labels.shape, dtype=bool)
    for name in ROAD_CLASSES:
        mask |= labels == CLASSES[name]
    # Schließen überbrückt die eingezeichneten Straßennamen.
    return closing(mask, disk(3))


# Wuchs der Landmaske in Pixeln (= Metern) für die Brückenprüfung. Am echten Mosaik
# gemessen: 3-5 px reichen nicht, weil Bordstein/Anti-Aliasing zwischen "land" und der
# angrenzenden Straße auf Karten mit echtem Rauschen mehr Lücke lassen als im
# künstlichen Test. Ist der Wuchs zu knapp, berührt eine ganz normale Küstenstraße
# nirgends "Land" und verkettet sich mit jeder anderen zu knappen Stelle zu einem
# Brückenstück, das (wie der ursprüngliche Fensteransatz) fast die ganze Karte
# überspannt. Ab 8 px reißt diese Kette am echten Mosaik; 10 px lässt zusätzlich Raum
# zum künstlichen Test, dessen 15 Pixel breite Wasserlücke sonst zu knapp anläge.
LAND_REACH_PX = 10


def road_polygons(labels: np.ndarray) -> list[tuple[Polygon, bool]]:
    """Straßenflächen, aufgeteilt in Stücke auf Land und Brückenstücke über Wasser.

    "Land" ist hier die Geländeklasse ohne die Straße selbst, um LAND_REACH_PX nach
    außen gewachsen. Ohne den Ausschluss der Straße wäre das zirkulär: land_polygons()
    zählt "nicht Wasser" als Land, und eine Straße ist nie "wasser" – ein Damm quer
    durchs Meer würde dann fast vollständig mit sich selbst überlappen und nie als
    Brücke auffallen. Der Wuchs sorgt dafür, dass eine Straße, die direkt an echtem
    Land beginnt, dort noch als "auf Land" erkannt wird, statt an ihrer eigenen Kante
    ins Leere zu greifen.
    """
    solid_ground = ndi.binary_dilation(
        (labels != CLASSES["wasser"]) & ~road_mask(labels), disk(LAND_REACH_PX)
    )
    land = unary_union(mask_polygons(solid_ground, min_area=50, simplify=1.0))
    out = []
    for poly in mask_polygons(road_mask(labels), min_area=20, simplify=0.8):
        for part, bridge in ((poly.intersection(land), False), (poly.difference(land), True)):
            for piece in getattr(part, "geoms", [part]):
                if piece.geom_type == "Polygon" and piece.area >= 20:
                    out.append((piece, bridge))
    return out


def centerlines(labels: np.ndarray) -> list[LineString]:
    """Skelett der Straßenflächen, nur zum Platzieren von Laternen und Palmen.

    Eine Kontur um das ein Pixel breite Skelett liefert jede Linie doppelt (sie umrundet
    das Band und läuft auf der anderen Seite zurück). Stattdessen wird das Skelett als
    Pixelgraph gelesen: Skelettpixel sind Knoten, 8er-Nachbarschaft ist eine Kante. Je
    Zusammenhangskomponente werden die Pfade zwischen Knoten mit Grad ungleich 2 verfolgt
    (Enden und Kreuzungen trennen die Äste); ein reiner Ring ohne solche Knoten wird
    einmal umlaufen.
    """
    graph = _skeleton_graph(skeletonize(road_mask(labels)))
    lines = []
    for path in _skeleton_paths(graph):
        line = LineString(_to_world(path)).simplify(1.5)
        if line.length >= 15:
            lines.append(line)
    return lines


def _skeleton_graph(skeleton: np.ndarray) -> nx.Graph:
    """Skelettpixel als Knoten, 8er-Nachbarschaft als Kanten."""
    pixels = set(zip(*np.nonzero(skeleton)))
    graph = nx.Graph()
    graph.add_nodes_from(pixels)
    for r, c in pixels:
        for dr in (-1, 0, 1):
            for dc in (-1, 0, 1):
                neighbour = (r + dr, c + dc)
                if (dr or dc) and neighbour in pixels:
                    graph.add_edge((r, c), neighbour)
    return graph


def _skeleton_paths(graph: nx.Graph) -> list[list[tuple[int, int]]]:
    """Pfade zwischen Enden/Kreuzungen (Grad ≠ 2) je Zusammenhangskomponente;
    ein reiner Ring aus lauter Grad-2-Knoten wird einmal umlaufen."""
    paths, seen = [], set()
    for component in nx.connected_components(graph):
        sub = graph.subgraph(component)
        branch_nodes = [n for n in sub if sub.degree(n) != 2]
        if not branch_nodes:
            start = next(iter(component))
            ring = list(nx.dfs_preorder_nodes(sub, start))
            ring.append(start)
            paths.append(ring)
            continue
        for start in branch_nodes:
            for first in sub.neighbors(start):
                if frozenset((start, first)) in seen:
                    continue
                path, prev, current = [start, first], start, first
                seen.add(frozenset((start, first)))
                while sub.degree(current) == 2:
                    nxt = next(n for n in sub.neighbors(current) if n != prev)
                    seen.add(frozenset((current, nxt)))
                    path.append(nxt)
                    prev, current = current, nxt
                paths.append(path)
    return paths