diff --git a/tools/extract/leonida/vectors.py b/tools/extract/leonida/vectors.py new file mode 100644 index 0000000..d7ab033 --- /dev/null +++ b/tools/extract/leonida/vectors.py @@ -0,0 +1,97 @@ +import numpy as np +from scipy import ndimage as ndi +from shapely.geometry import LineString, Polygon +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: + 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)) + + +BRIDGE_MARGIN_PX = 3 # Rand um die Straßenfläche; die Fläche selbst ist nie "wasser". + + +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) + 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))) + 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)) + 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) + return lines diff --git a/tools/extract/tests/test_vectors.py b/tools/extract/tests/test_vectors.py new file mode 100644 index 0000000..b6c2da0 --- /dev/null +++ b/tools/extract/tests/test_vectors.py @@ -0,0 +1,33 @@ +import numpy as np +from leonida.classify import CLASSES +from leonida.vectors import land_polygons, mask_polygons, road_polygons + +def test_square_becomes_one_polygon_with_right_area(): + mask = np.zeros((100, 100), dtype=bool) + mask[20:60, 30:80] = True # 40 x 50 Pixel = 2000 m² + polys = mask_polygons(mask, min_area=10, simplify=1.0) + assert len(polys) == 1 + assert abs(polys[0].area - 2000) < 60 # Vereinfachung darf etwas abweichen + +def test_island_keeps_its_lagoon(): + labels = np.full((80, 80), CLASSES["wasser"], dtype=np.uint8) + labels[10:70, 10:70] = CLASSES["land"] + labels[30:50, 30:50] = CLASSES["wasser"] # Lagune + polys = land_polygons(labels) + assert len(polys) == 1 + assert len(polys[0].interiors) == 1 + +def test_road_over_water_is_flagged_as_bridge(): + labels = np.full((60, 60), CLASSES["wasser"], dtype=np.uint8) + labels[0:10, :] = CLASSES["land"] + labels[25:31, :] = CLASSES["strasse"] # Damm quer durchs Wasser + roads = road_polygons(labels) + assert len(roads) == 1 + assert roads[0][1] is True + +def test_crossing_roads_form_a_single_polygon(): + labels = np.full((60, 60), CLASSES["land"], dtype=np.uint8) + labels[28:32, :] = CLASSES["strasse"] + labels[:, 28:32] = CLASSES["strasse"] + roads = road_polygons(labels) + assert len(roads) == 1 # Kreuzung, kein Zerfall in vier Arme