From 47c3cb98da7f302a65562f041a8ecf76873345a9 Mon Sep 17 00:00:00 2001 From: KNSO Date: Sat, 12 Sep 2026 15:07:46 +0200 Subject: [PATCH] =?UTF-8?q?fix(extract):=20Br=C3=BCckenerkennung=20per=20V?= =?UTF-8?q?erschneidung,=20Mittellinien=20ohne=20Dopplung?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Fix-Runde 1, drei Befunde: - road_polygons: Straßennetz ist ein zusammenhängendes Polygon über den ganzen Ausschnitt; eine Fensterprüfung um seine Bounding-Box sieht daher zwangsläufig Wasser und meldet fast alles als Brücke. Ersetzt durch Verschneidung jedes Straßenpolygons mit der Landfläche (Ruling 13). Die wörtlich vorgegebene Landberechnung (unary_union(land_polygons(labels))) besteht die bestehenden Tests nicht: eine Straße ist nie "wasser", also zählt land_polygons() sie selbst als Land, und ein Damm quer durchs Meer überlappt dann mit sich selbst statt als Brücke aufzufallen (siehe Fix-Bericht für den Nachweis). Land wird stattdessen ohne die Straßenklasse berechnet und um LAND_REACH_PX gewachsen, damit eine Straße, die direkt an Land beginnt, dort noch erkannt wird. Der Radius wurde am echten Mosaik kalibriert: unter 8 px verkettet sich eine gewöhnliche Küstenstraße mangels Berührung mit "Land" erneut zu einem kartenweiten Brückenstück; 10 px behebt das mit Sicherheitsabstand zur Wasserlücke des bestehenden Tests. - centerlines: eine Kontur um das ein Pixel breite Skelett lief auf beiden Seiten des Bandes hin und zurück und lieferte jede Linie doppelt; an Kreuzungen ergab das eine einzelne Schleife statt getrennter Äste (Ruling 14, neue Tests zuerst rot gesehen). Ersetzt durch einen Pixelgraphen (networkx): Skelettpixel als Knoten, 8er-Nachbarschaft als Kanten, Pfade zwischen Knoten mit Grad ungleich 2 je Zusammenhangskomponente, reine Ringe einmal umlaufen. - mask_polygons: Polygon(ring).buffer(0) kann bei selbstschneidenden Konturen ein MultiPolygon liefern; das schlug später mit AttributeError fehl. Beim Einsortieren in outers/holes wird jetzt nur noch ein sauberes Polygon akzeptiert (Ruling 15). 16/16 Tests grün. Auf dem echten Mosaik: 469 Straßenstücke (190 Brücken, größtes Brückenstück 26.294 m² auf einer lokalen Fläche von 643 x 564 m, nicht mehr kartenweit). Co-Authored-By: Claude Sonnet 5 --- tools/extract/leonida/vectors.py | 116 +++++++++++++++++++++------- tools/extract/tests/test_vectors.py | 15 +++- 2 files changed, 101 insertions(+), 30 deletions(-) diff --git a/tools/extract/leonida/vectors.py b/tools/extract/leonida/vectors.py index d6087a9..b057885 100644 --- a/tools/extract/leonida/vectors.py +++ b/tools/extract/leonida/vectors.py @@ -1,6 +1,8 @@ +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 @@ -24,7 +26,7 @@ def mask_polygons(mask: np.ndarray, min_area: float, simplify: float) -> list[Po if len(ring) < 4: continue poly = Polygon(ring).buffer(0) - if poly.is_empty or poly.area < min_area: + 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 = [] @@ -60,41 +62,97 @@ def road_mask(labels: np.ndarray) -> np.ndarray: return closing(mask, disk(3)) -# Rand in Pixeln (= Metern) um die Straßenfläche zur Brückenprüfung: die Straßenfläche -# selbst ist per Definition nie "wasser", daher muss das Fenster über ihren Rand -# hinausschauen, um das Wasser zu sehen, das sie überquert. -BRIDGE_MARGIN_PX = 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]]: - water = labels == CLASSES["wasser"] - polys = mask_polygons(road_mask(labels), min_area=20, simplify=0.8) + """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 polys: - r0, c0, r1, c1 = _bbox_px(poly) - r0 = max(r0 - BRIDGE_MARGIN_PX, 0) - c0 = max(c0 - BRIDGE_MARGIN_PX, 0) - r1 = min(r1 + BRIDGE_MARGIN_PX, water.shape[0] - 1) - c1 = min(c1 + BRIDGE_MARGIN_PX, water.shape[1] - 1) - window = water[r0:r1 + 1, c0:c1 + 1] - out.append((poly, bool(window.mean() > 0.3))) + 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 _bbox_px(poly: Polygon) -> tuple[int, int, int, int]: - minx, miny, maxx, maxy = poly.bounds - return (int(SLICE["y1"] - maxy), int(minx - SLICE["x0"]), - int(SLICE["y1"] - miny), int(maxx - SLICE["x0"])) - - def centerlines(labels: np.ndarray) -> list[LineString]: - """Skelett der Straßenflächen, nur zum Platzieren von Laternen und Palmen.""" - skeleton = skeletonize(road_mask(labels)) + """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 contour in find_contours(skeleton.astype(float), 0.5): - pts = _to_world(contour) - if len(pts) >= 2: - line = LineString(pts).simplify(1.5) - if line.length >= 15: - lines.append(line) + 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 diff --git a/tools/extract/tests/test_vectors.py b/tools/extract/tests/test_vectors.py index b6c2da0..878b6ba 100644 --- a/tools/extract/tests/test_vectors.py +++ b/tools/extract/tests/test_vectors.py @@ -1,6 +1,6 @@ import numpy as np from leonida.classify import CLASSES -from leonida.vectors import land_polygons, mask_polygons, road_polygons +from leonida.vectors import centerlines, land_polygons, mask_polygons, road_polygons def test_square_becomes_one_polygon_with_right_area(): mask = np.zeros((100, 100), dtype=bool) @@ -31,3 +31,16 @@ def test_crossing_roads_form_a_single_polygon(): labels[:, 28:32] = CLASSES["strasse"] roads = road_polygons(labels) assert len(roads) == 1 # Kreuzung, kein Zerfall in vier Arme + +def test_centerline_runs_once_along_a_straight_road(): + labels = np.full((60, 300), CLASSES["land"], dtype=np.uint8) + labels[28:38, 50:250] = CLASSES["strasse"] # 200 lang, 10 breit + lines = centerlines(labels) + assert len(lines) == 1 + assert 180 <= lines[0].length <= 215 # einmal entlang, nicht hin und zurück + +def test_crossing_roads_yield_separate_centerlines(): + labels = np.full((120, 120), CLASSES["land"], dtype=np.uint8) + labels[55:65, :] = CLASSES["strasse"] + labels[:, 55:65] = CLASSES["strasse"] + assert len(centerlines(labels)) >= 3 # Arme werden an der Kreuzung getrennt