gta6/tools/extract/leonida/roads.py

"""Straßennetz als Graph statt als nachgezeichnete Kontur.

Der Entwurf (docs/superpowers/specs/2026-09-16-strassen-design.md) misst nach, dass die
Mittellinien im Median nur um den Faktor 1,0104 von der Luftlinie abweichen: Die Straßen sind
längst gerade, das Zackige sitzt allein in den Rändern von road_polygons(). Wer die Ränder aus
Mittellinie und Breite neu baut, statt sie zu glätten, bekommt gerade Straßen fast geschenkt.

centerlines() in vectors.py taugt dafür nicht: Sein Filter length >= 15 wirft an jeder Kreuzung
die kurzen Stummel weg, und von 3.038 Linienenden treffen sich dadurch nur 81 zu Kreuzungen.
Hier bleibt das Skelett deshalb ungefiltert, und die Kreuzungen entstehen aus den Pixeln mit
Grad >= 3 statt aus zufällig aufeinandertreffenden Linienenden.
"""
from __future__ import annotations

from collections import Counter
from dataclasses import dataclass, field

import networkx as nx
import numpy as np
from scipy import ndimage as ndi
from shapely.geometry import LineString, Point, Polygon
from shapely.ops import substring
from skimage.morphology import disk

from .coords import SLICE
from .vectors import skeleton_graph, to_world

# Radius, in dem Kreuzungspixel zu einem Knoten verschmelzen. Am echten Mosaik gemessen: Die
# 5.970 Pixel mit Grad >= 3 fallen damit auf 862 Knoten zusammen. Kleiner, und eine breite
# Kreuzung zerfällt in mehrere Knoten, weil das Skelett dort sternförmig ausfranst.
JUNCTION_RADIUS = 8


@dataclass(frozen=True)
class Node:
    point: tuple[float, float]


@dataclass
class Edge:
    a: int
    b: int
    points: list[tuple[float, float]]
    width: float = 0.0
    klass: str = ""


@dataclass
class Graph:
    nodes: list[Node] = field(default_factory=list)
    edges: list[Edge] = field(default_factory=list)

    def incident(self, node: int) -> list[int]:
        return [k for k, e in enumerate(self.edges) if e.a == node or e.b == node]

    def degree(self, node: int) -> int:
        # Eine Schleife (a == b) zählt zweimal, wie in jedem Graphen.
        return sum((e.a == node) + (e.b == node) for e in self.edges)


def graph_from_skeleton(skel: np.ndarray) -> Graph:
    """Knoten aus den Kreuzungspixeln, Kanten aus dem, was zwischen ihnen übrig bleibt."""
    pixels = skeleton_graph(skel)
    clusters, points = _junction_clusters(pixels, skel.shape)
    graph = Graph(nodes=[Node(p) for p in points])
    rest = pixels.copy()
    rest.remove_nodes_from([p for p in pixels if clusters[p] > 0])
    for component in nx.connected_components(rest):
        chain = _chain(rest.subgraph(component))
        if not chain:
            continue
        a = _attach(graph, pixels, clusters, chain[0])
        b = _attach(graph, pixels, clusters, chain[-1])
        line = to_world(chain)
        if graph.nodes[a].point != line[0]:
            line.insert(0, graph.nodes[a].point)
        if graph.nodes[b].point != line[-1]:
            line.append(graph.nodes[b].point)
        graph.edges.append(Edge(a=a, b=b, points=line))
    return graph


def _junction_clusters(pixels: nx.Graph, shape) -> tuple[np.ndarray, list[tuple[float, float]]]:
    """Kreuzungspixel, um JUNCTION_RADIUS gewachsen, in Zusammenhangskomponenten.

    Der Knotenpunkt ist der Schwerpunkt der *ungewachsenen* Pixel: Das Wachsen dient nur dazu,
    benachbarte Kreuzungen zusammenzufassen, und würde den Punkt sonst nach außen ziehen.
    """
    junction = np.zeros(shape, dtype=bool)
    for pixel in pixels:
        if pixels.degree(pixel) >= 3:
            junction[pixel] = True
    labelled, count = ndi.label(ndi.binary_dilation(junction, disk(JUNCTION_RADIUS)))
    points = []
    for k in range(1, count + 1):
        rows, cols = np.nonzero(junction & (labelled == k))
        points.append(to_world([(rows.mean(), cols.mean())])[0])
    # labelled traegt die Knotennummer auf jedem *gewachsenen* Pixel. Genau das braucht _attach():
    # Das Ende einer Kette liegt neben einem gewachsenen Pixel, nicht neben einem Kreuzungspixel.
    return labelled, points


def _chain(sub: nx.Graph) -> list[tuple[int, int]]:
    """Die Pixel einer Kette in Laufrichtung. Ohne Kreuzungspixel ist jede Komponente eine."""
    ends = [n for n in sub if sub.degree(n) <= 1]
    if not ends:
        return []            # geschlossener Ring ohne Anschluss: keine Straße, die irgendwo hinführt
    if len(ends) == 1:
        return [ends[0]]
    return nx.shortest_path(sub, ends[0], ends[-1])


def _attach(graph: Graph, pixels: nx.Graph, clusters: np.ndarray, pixel) -> int:
    """Knoten-ID am Ende einer Kette: der berührte Kreuzungsknoten oder ein neues Sackgassenende."""
    for neighbour in pixels.neighbors(pixel):
        if clusters[neighbour] > 0:
            return int(clusters[neighbour]) - 1
    graph.nodes.append(Node(to_world([pixel])[0]))
    return len(graph.nodes) - 1


# Sackgassen unter dieser Länge sind Skelettfransen an breiten Flächen, keine Straßen.
SPUR_MIN = 10.0
# Zwei Kreuzungen näher als das sind eine: das Skelett franst an Parkplätzen sternförmig aus.
CONTRACT_MAX = 5.0
SIMPLIFY = 1.0
# Die ersten und letzten Meter einer Kante liegen im Kreuzungsbereich; dort misst die
# Distanztransformation die Breite der Kreuzung, nicht die der Straße.
WIDTH_MARGIN = 8.0
WIDTH_STEP = 5.0
CLASS_WIDTH = {"gasse": 6.0, "wohnstrasse": 12.0, "hauptstrasse": 18.0, "boulevard": 24.0}


def road_graph(labels: np.ndarray, distance: np.ndarray) -> Graph:
    """Der Aufruf für die Pipeline: Skelett, Graph, säubern, Breite und Klasse."""
    from skimage.morphology import skeletonize

    from .vectors import road_mask

    return measure_widths(clean(graph_from_skeleton(skeletonize(road_mask(labels)))), distance)


def clean(graph: Graph) -> Graph:
    """Fransen weg, Zwischenknoten auflösen, zu kurze Verbindungen zusammenziehen, glätten."""
    out = Graph(nodes=list(graph.nodes), edges=list(graph.edges))
    _drop_spurs(out)
    _merge_chains(out)
    _contract_short(out)
    _merge_chains(out)
    for edge in out.edges:
        simple = LineString(edge.points).simplify(SIMPLIFY)
        edge.points = [(round(x, 2), round(y, 2)) for x, y in simple.coords]
    return _renumber(out)


def _length(edge: Edge) -> float:
    # Ein Rohast ohne Kreuzungsanschluss liefert graph_from_skeleton() als Kante mit nur einem
    # Stützpunkt (Anfang und Ende landen auf demselben neuen Sackgassenknoten) - LineString
    # verlangt 0 oder mehr als einen Punkt. Geometrisch ist das kein Weg, also Länge 0: Beide
    # Enden bleiben Grad 1, und die Fransenregel entfernt ihn ganz regulär im nächsten Durchgang.
    return LineString(edge.points).length if len(edge.points) > 1 else 0.0


def _degrees(graph: Graph) -> Counter:
    grade = Counter()
    for edge in graph.edges:
        grade[edge.a] += 1
        grade[edge.b] += 1
    return grade


def _drop_spurs(graph: Graph) -> None:
    """Fransen weg, in Durchgängen statt Fund für Fund.

    Der Grad wird einmal gezählt und beim Entfernen fortgeschrieben, statt ihn je Kante neu über
    alle Kanten zu suchen, und ein Durchgang räumt alle Fransen auf einmal weg, statt nach jedem
    Fund von vorn zu beginnen. Beides zusammen macht aus einer quadratischen eine lineare
    Rechnung je Durchgang — das echte Skelett hat rund 12.000 Rohäste, und die naive Fassung
    rechnete daran Stunden. Die Zahl der Durchgänge ist die Tiefe der Fransenketten, also klein.
    """
    grade = _degrees(graph)
    while True:
        bleibt, entfernt = [], 0
        for edge in graph.edges:
            if (grade[edge.a] == 1 or grade[edge.b] == 1) and _length(edge) < SPUR_MIN:
                grade[edge.a] -= 1
                grade[edge.b] -= 1
                entfernt += 1
            else:
                bleibt.append(edge)
        graph.edges = bleibt
        if not entfernt:
            return


def _incident_index(graph: Graph) -> dict[int, list[int]]:
    index: dict[int, list[int]] = {}
    for k, edge in enumerate(graph.edges):
        index.setdefault(edge.a, []).append(k)
        if edge.b != edge.a:
            index.setdefault(edge.b, []).append(k)
    return index


def _merge_chains(graph: Graph) -> None:
    """Knoten mit genau zwei Kanten auflösen. Der Index entsteht je Verschmelzung einmal neu —
    das ist vertretbar, weil es nach dem Fransenschnitt nur noch wenige solcher Knoten gibt."""
    changed = True
    while changed:
        changed = False
        for node, incident in _incident_index(graph).items():
            if len(incident) != 2:
                continue
            first, second = (graph.edges[i] for i in incident)
            if first.a == first.b or second.a == second.b:
                continue          # eine Schleife (z.B. ein kleiner Ring am Skelett) hat kein
                                  # zweites Ende - mit der Nachbarkante lässt sie sich nicht verketten
            points = _join(first, second, graph.nodes[node].point)
            if points is None:
                continue
            ends = [n for e in (first, second) for n in (e.a, e.b) if n != node]
            for i in sorted(incident, reverse=True):
                del graph.edges[i]
            graph.edges.append(Edge(a=ends[0], b=ends[1], points=points))
            changed = True
            break


def _join(first: Edge, second: Edge, at: tuple[float, float]) -> list[tuple[float, float]] | None:
    """Zwei Kanten an einem gemeinsamen Punkt zu einer Punktfolge verketten."""
    a = first.points if first.points[-1] == at else list(reversed(first.points))
    b = second.points if second.points[0] == at else list(reversed(second.points))
    if a[-1] != at or b[0] != at:
        return None          # die Enden passen nicht zusammen: nichts anfassen
    return a + b[1:]


def _contract_short(graph: Graph) -> None:
    changed = True
    while changed:
        changed = False
        grade = _degrees(graph)
        for k, edge in enumerate(graph.edges):
            if edge.a == edge.b or _length(edge) >= CONTRACT_MAX:
                continue
            if grade[edge.a] < 3 or grade[edge.b] < 3:
                continue
            ax, ay = graph.nodes[edge.a].point
            bx, by = graph.nodes[edge.b].point
            graph.nodes[edge.a] = Node(((ax + bx) / 2, (ay + by) / 2))
            del graph.edges[k]
            for other in graph.edges:
                if other.a == edge.b:
                    other.a = edge.a
                if other.b == edge.b:
                    other.b = edge.a
                # Auch die Kanten, die schon vorher am überlebenden Knoten hingen, müssen ihm
                # folgen: Er ist gerade auf den Mittelpunkt zwischen beiden Kreuzungen gewandert.
                if other.a == edge.a:
                    other.points = [graph.nodes[edge.a].point, *other.points[1:]]
                if other.b == edge.a:
                    other.points = [*other.points[:-1], graph.nodes[edge.a].point]
            changed = True
            break


def _renumber(graph: Graph) -> Graph:
    """Knoten ohne Kante fallen weg; die Kanten zeigen danach auf die neuen Nummern."""
    used = sorted({n for e in graph.edges for n in (e.a, e.b)})
    index = {old: new for new, old in enumerate(used)}
    return Graph(nodes=[graph.nodes[i] for i in used],
                 edges=[Edge(index[e.a], index[e.b], e.points, e.width, e.klass) for e in graph.edges])


def measure_widths(graph: Graph, distance: np.ndarray) -> Graph:
    """Breite je Kante aus der Distanztransformation, dann die nächstliegende Regelbreite."""
    for edge in graph.edges:
        line = LineString(edge.points)
        edge.width = _median_width(line, distance)
        edge.klass = min(CLASS_WIDTH, key=lambda name: abs(CLASS_WIDTH[name] - edge.width))
    return graph


# Kürzere Reste nach einem Schnitt sind Stummel, die nur eine Kreuzung vortäuschen würden.
MIN_EDGE = 5.0


def cut_edges(graph: Graph, cut: Polygon) -> Graph:
    """Alles, was in der Fläche liegt, aus den Kanten herausschneiden.

    Gegenstück zum Schneiden der Straßenflächen: Wo Flamingo Park die alte Gasse entfernt, darf
    im Graphen an derselben Stelle keine Kante mehr stehen - sonst führte dort weiterhin eine
    Fahrbahn durch den Park, nur ohne die Straßenfläche, die sie einmal begründet hat.
    """
    out = Graph(nodes=list(graph.nodes), edges=[])
    for edge in graph.edges:
        line = LineString(edge.points)
        if not line.intersects(cut):
            out.edges.append(edge)
            continue
        rest = line.difference(cut)
        for piece in getattr(rest, "geoms", [rest]):
            if piece.geom_type != "LineString" or piece.length < MIN_EDGE:
                continue
            points = [(round(x, 2), round(y, 2)) for x, y in piece.coords]
            out.edges.append(Edge(a=_cut_end(out, points[0], edge), b=_cut_end(out, points[-1], edge),
                                  points=points, width=edge.width, klass=edge.klass))
    return _renumber(out)


def _cut_end(graph: Graph, point: tuple[float, float], edge: Edge) -> int:
    """Ein Ende eines Schnittstücks: der ursprüngliche Knoten, wenn das Stück dort endet, sonst ein
    neues Sackgassenende. Nicht attach(): dessen 8-m-Umkreis fände bei einem Reststück unter 8 m für
    beide Enden denselben Knoten und machte aus dem Stück eine Schleife."""
    for node in (edge.a, edge.b):
        px, py = graph.nodes[node].point
        if abs(px - point[0]) < 0.01 and abs(py - point[1]) < 0.01:
            return node
    graph.nodes.append(Node((float(point[0]), float(point[1]))))
    return len(graph.nodes) - 1


def attach(graph: Graph, point: tuple[float, float], radius: float = JUNCTION_RADIUS) -> int:
    """Knoten-ID an dieser Stelle: der nächste innerhalb von radius, sonst ein neuer.

    Damit hängt sich eine von Hand ergänzte Straße (die Jefferson Avenue) an die Querstraße an,
    statt als freies Ende danebenzuliegen und ohne Kreuzung zu enden.
    """
    best, best_distance = None, radius
    for k, node in enumerate(graph.nodes):
        distance = ((node.point[0] - point[0]) ** 2 + (node.point[1] - point[1]) ** 2) ** 0.5
        if distance <= best_distance:
            best, best_distance = k, distance
    if best is not None:
        return best
    graph.nodes.append(Node((float(point[0]), float(point[1]))))
    return len(graph.nodes) - 1


# Kürzestes Teilstück, das insert_edge() beim Aufbrechen einer gekreuzten Kante stehen lässt;
# kürzer wäre praktisch derselbe Punkt und keine Fahrbahn. Dieselbe Grenze wie trimmed_line() in
# section.py für "die Kante liegt ganz in der Kreuzung".
INSERT_MIN = 1.0


def insert_edge(graph: Graph, points: list[tuple[float, float]], width: float, klass: str) -> Graph:
    """Eine von Hand ergänzte Straße so ins Netz einfügen, dass sie an ihren Kreuzungen Knoten bildet.

    attach() allein reicht nur, wenn ein Ende der neuen Straße neben einem vorhandenen Knoten
    liegt. Die Jefferson Avenue lief am echten Netz mittendrin über die 11th Street hinweg, ohne
    eines ihrer beiden Enden dort zu haben: attach() hätte die neue Fahrbahn einfach über die alte
    gelegt, ohne dass an der Stelle eine Kreuzungsfläche, ein Eckradius, ein Zebrastreifen oder
    eine Haltelinie entsteht.
    """
    out = Graph(nodes=list(graph.nodes), edges=[])
    new_line = LineString(points)

    # Jede vorhandene Kante, die die neue Linie kreuzt, an den Kreuzungspunkten aufbrechen.
    # attach() entscheidet je Punkt, ob dort ein neuer Knoten entsteht oder ein vorhandener in
    # der Nähe genügt (Regel b) - dieselbe Funktion, die gleich auch die Enden der neuen Straße
    # anhängt, damit an einer Kreuzung nicht zwei Knoten wenige Meter auseinander liegen.
    crossings: list[tuple[float, float]] = []
    for edge in graph.edges:
        line = LineString(edge.points)
        found = _crossing_points(new_line, line)
        if not found:
            out.edges.append(edge)
            continue
        crossings += found
        stops = sorted({0.0, line.length} | {line.project(Point(p)) for p in found})
        for start, end in zip(stops, stops[1:]):
            _add_piece(out, line, start, end, edge.width, edge.klass)

    # Dieselben Kreuzungspunkte (jetzt echte Knoten) plus die beiden Enden der neuen Straße.
    stops = sorted({0.0, new_line.length} | {new_line.project(Point(p)) for p in crossings})
    for start, end in zip(stops, stops[1:]):
        _add_piece(out, new_line, start, end, width, klass)
    return _renumber(out)


def _crossing_points(a: LineString, b: LineString) -> list[tuple[float, float]]:
    """Die Punkte, an denen sich zwei Linien berühren - deckungsgleiche Teilstücke (kein Punkt,
    sondern selbst eine Linie) sind ein Datenfehler und werden hier ignoriert, nicht aufgebrochen."""
    hit = a.intersection(b)
    if hit.geom_type == "Point":
        return [(hit.x, hit.y)]
    if hit.geom_type == "MultiPoint":
        return [(p.x, p.y) for p in hit.geoms]
    return []


def _add_piece(out: Graph, line: LineString, start: float, end: float, width: float, klass: str) -> None:
    """Das Teilstück von `line` zwischen zwei Stationen anhängen, wenn es lang genug ist und seine
    Enden nicht auf denselben Knoten fallen (Regel d) - attach() sucht dafür je Ende in `out`."""
    if end - start < INSERT_MIN:
        return
    points = [(round(x, 2), round(y, 2)) for x, y in substring(line, start, end).coords]
    a, b = attach(out, points[0]), attach(out, points[-1])
    if a == b:
        return
    out.edges.append(Edge(a=a, b=b, points=points, width=width, klass=klass))


def _median_width(line: LineString, distance: np.ndarray) -> float:
    start, end = WIDTH_MARGIN, line.length - WIDTH_MARGIN
    if end <= start:
        start, end = 0.0, line.length      # zu kurz für die Randregel: dann eben ganz messen
    steps = max(1, int((end - start) // WIDTH_STEP))
    samples = []
    for i in range(steps + 1):
        point = line.interpolate(start + (end - start) * i / steps)
        row = int(round(SLICE["y1"] - point.y))
        col = int(round(point.x - SLICE["x0"]))
        if 0 <= row < distance.shape[0] and 0 <= col < distance.shape[1]:
            samples.append(distance[row, col] * 2)
    return float(np.median(samples)) if samples else 0.0