feat(extract): Graph säubern, Breite messen, Klasse bestimmen

Sackgassen unter 10 m sind Skelettfransen und fallen weg, Knoten mit nur
noch zwei Kanten lösen sich auf, Kreuzungen näher als 5 m werden eine. Die
Breite kommt aus der Distanztransformation ohne die ersten und letzten 8 m:
dort misst sie die Kreuzung, nicht die Straße.

Der Gegencheck am echten Mosaik deckte zwei Fälle auf, die kein
synthetischer Test zeigt: Ein Rohast ohne Kreuzungsanschluss liefert eine
Kante mit nur einem Stützpunkt (LineString verlangt 0 oder mehr als einen
Punkt), und ein Knoten mit Schleife plus Verbindungskante hat kein zweites
Ende zum Verketten. Beides fängt clean() jetzt ab, mit je einem
Regressionstest.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>

tools/extract/leonida/roads.py

@@ -12,13 +12,16 @@ 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
 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
@@ -113,3 +116,176 @@ def _attach(graph: Graph, pixels: nx.Graph, clusters: np.ndarray, pixel) -> int:
             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
+                    other.points = [graph.nodes[edge.a].point, *other.points[1:]]
+                if other.b == edge.b:
+                    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
+
+
+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