diff --git a/tools/extract/leonida/buildings.py b/tools/extract/leonida/buildings.py new file mode 100644 index 0000000..c6c9b1e --- /dev/null +++ b/tools/extract/leonida/buildings.py @@ -0,0 +1,40 @@ +import warnings + +import numpy as np +from shapely.geometry import Polygon + +from .classify import CLASSES +from .vectors import mask_polygons + +MIN_AREA = 12.0 +RECT_COVERAGE = 0.85 + + +def regularize(poly: Polygon) -> Polygon | None: + """Rechteck, wenn das Polygon sein gedrehtes Hüllrechteck gut ausfüllt; sonst vereinfacht.""" + poly = poly.buffer(0) + if poly.geom_type != "Polygon" or poly.area < MIN_AREA: + return None + with warnings.catch_warnings(): + # Shapelys oriented_envelope (minimum_rotated_rectangle) wertet einen Winkelterm + # aus, der für bereits achsenparallele/rechtwinklige Polygone exakt 0/0 wird - + # GEOS meldet das als RuntimeWarning, obwohl das zurückgegebene Rechteck korrekt + # bleibt (geprüft: Fläche stimmt, kein NaN). Nur diese eine Property-Abfrage ist + # eingeschränkt, damit ein echter Warnhinweis an anderer Stelle sichtbar bliebe. + warnings.filterwarnings("ignore", message=r".*oriented_envelope.*", category=RuntimeWarning) + rect = poly.minimum_rotated_rectangle + if rect.area > 0 and poly.area / rect.area >= RECT_COVERAGE: + return rect + simple = poly.simplify(1.0) + return simple if simple.area >= MIN_AREA else None + + +def building_polygons(labels: np.ndarray) -> list[tuple[Polygon, bool]]: + out = [] + for name, speculative in (("gebaeude", False), ("spekulativ", True)): + mask = labels == CLASSES[name] + for poly in mask_polygons(mask, min_area=MIN_AREA, simplify=0.5): + reg = regularize(poly) + if reg is not None: + out.append((reg, speculative)) + return out diff --git a/tools/extract/leonida/heights.py b/tools/extract/leonida/heights.py new file mode 100644 index 0000000..86732d7 --- /dev/null +++ b/tools/extract/leonida/heights.py @@ -0,0 +1,102 @@ +import hashlib +import math + +from shapely.geometry import Point, Polygon + +DISTRICT_RANGES = { + "Downtown": (20.0, 60.0), + "Ocean Beach": (8.0, 25.0), + "default": (6.0, 15.0), +} +LANDMARK_SNAP = 20.0 # Meter: so weit darf ein Wahrzeichenpunkt vom Grundriss entfernt liegen + + +def _seeded_unit(seed: int, key: str) -> float: + digest = hashlib.sha256(f"{seed}:{key}".encode()).digest() + return int.from_bytes(digest[:8], "big") / 2 ** 64 + + +def assign_heights(polys: list[Polygon], landmarks: list[dict], + districts: list[dict], seed: int) -> list[dict]: + out = [] + for poly in polys: + centroid = poly.centroid + key = f"{centroid.x:.1f},{centroid.y:.1f}" + hit = None + for lm in landmarks: + p = Point(lm["xyz"][0], lm["xyz"][1]) + if poly.contains(p) or poly.distance(p) <= LANDMARK_SNAP: + if hit is None or lm["height"] > hit["height"]: + hit = lm + if hit is not None: + out.append({"height": float(hit["height"]), "source": "skyline"}) + continue + name = "default" + best = math.inf + for d in districts: + dist = math.dist((centroid.x, centroid.y), (d["x"], d["y"])) + if dist < best: + best, name = dist, d["name"] if d["name"] in DISTRICT_RANGES else "default" + lo, hi = DISTRICT_RANGES[name] + # Kleine Grundrisse bleiben niedriger: Fläche skaliert den oberen Rand. + area_factor = min(1.0, poly.area / 400.0) + height = lo + (hi - lo) * _seeded_unit(seed, key) * (0.4 + 0.6 * area_factor) + out.append({"height": round(height, 1), "source": "district"}) + return out + + +PALM_SPACING = 12.0 # Meter entlang der Mittellinie +PALM_BEACH_RANGE = 40.0 # Meter Abstand zur nächsten Sandfläche +PALM_CLEARANCE = 6.0 # Meter Abstand zu Gebäuden +PALM_KERB = 1.5 # Meter Abstand zum Fahrbahnrand +# Versatzkandidaten von innen nach außen; genommen wird der erste, der Fahrbahn und +# Gebäude freihält. Ein fester Wert geht nicht: die Fahrbahn ist an den Mittellinien- +# punkten p50 5,4 m, p90 11,7 m, p95 14,0 m halbbreit (am echten Mosaik gemessen). +PALM_VERGE = (6.0, 8.0, 10.0, 12.0, 14.0, 16.0) + + +def _verge_point(line, p, side: str, roads, buildings): + """Erster Versatz von innen nach außen, der Fahrbahn und Gebäude freihält.""" + for dist in PALM_VERGE: + off = line.parallel_offset(dist, side) + if off.is_empty: + continue + q = off.interpolate(off.project(p)) + if roads is not None and q.distance(roads) < PALM_KERB: + continue + if buildings is not None and q.distance(buildings) < PALM_CLEARANCE: + continue + return [(round(q.x, 2), round(q.y, 2))] + return [] + + +def palm_points(sand_polys, park_polys, road_lines, building_polys, road_polys, seed: int): + """Palmen entlang strandnaher Mittellinien und verstreut in Parks.""" + from shapely.ops import unary_union + sand = unary_union(sand_polys) if sand_polys else None + roads = unary_union(list(road_polys)) if road_polys else None + buildings = unary_union(list(building_polys)) if building_polys else None + blocked = unary_union([g for g in (roads, buildings) if g is not None]) \ + if (roads is not None or buildings is not None) else None + points = [] + for line in road_lines: + if sand is None or line.distance(sand) > PALM_BEACH_RANGE: + continue + # Start bei halber Teilung: die Mittellinien sind an Kreuzungen getrennt, ein + # Start bei 0 setzt je Ast eine Palme in denselben Punkt. + pos = min(PALM_SPACING / 2, line.length / 2) + while pos <= line.length: + p = line.interpolate(pos) + for side in ("left", "right"): + points.extend(_verge_point(line, p, side, roads, buildings)) + pos += PALM_SPACING + for park in park_polys: + minx, miny, maxx, maxy = park.bounds + n = int(park.area // 250) + for k in range(n): + u = _seeded_unit(seed, f"park:{minx:.0f}:{k}:x") + v = _seeded_unit(seed, f"park:{miny:.0f}:{k}:y") + q = Point(minx + u * (maxx - minx), miny + v * (maxy - miny)) + if park.contains(q) and (blocked is None or q.distance(blocked) > PALM_CLEARANCE): + points.append((round(q.x, 2), round(q.y, 2))) + return points diff --git a/tools/extract/pyproject.toml b/tools/extract/pyproject.toml index a03e290..720e769 100644 --- a/tools/extract/pyproject.toml +++ b/tools/extract/pyproject.toml @@ -13,3 +13,11 @@ dev = ["pytest>=8.3"] [tool.pytest.ini_options] testpaths = ["tests"] pythonpath = ["."] +filterwarnings = [ + # Shapelys oriented_envelope (minimum_rotated_rectangle) meldet fuer bereits + # achsenparallele/rechtwinklige Polygone eine RuntimeWarning aus GEOS (Winkelterm + # wird 0/0), obwohl das zurueckgegebene Rechteck nachweislich korrekt bleibt. Nur + # diese eine, an Text und Kategorie festgemachte Warnung wird gefiltert; jede andere + # Warnung (auch jede echte DeprecationWarning) bleibt sichtbar. + "ignore:.*oriented_envelope.*:RuntimeWarning", +] diff --git a/tools/extract/tests/test_buildings.py b/tools/extract/tests/test_buildings.py new file mode 100644 index 0000000..dd877c0 --- /dev/null +++ b/tools/extract/tests/test_buildings.py @@ -0,0 +1,24 @@ +import numpy as np +from shapely.geometry import Polygon +from leonida.buildings import building_polygons, regularize +from leonida.classify import CLASSES + +def test_rough_rectangle_snaps_to_rectangle(): + poly = Polygon([(0, 0), (10.4, 0.3), (10.1, 20.2), (-0.2, 19.8)]) + out = regularize(poly) + assert len(out.exterior.coords) == 5 # geschlossenes Rechteck + assert abs(out.area - 205) < 25 + +def test_l_shape_is_kept_not_boxed(): + poly = Polygon([(0, 0), (30, 0), (30, 10), (10, 10), (10, 30), (0, 30)]) + out = regularize(poly) + assert out.area < 0.8 * out.minimum_rotated_rectangle.area + +def test_tiny_polygon_is_dropped(): + assert regularize(Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])) is None + +def test_two_buildings_are_separated(): + labels = np.full((60, 60), CLASSES["stadt"], dtype=np.uint8) + labels[10:25, 10:25] = CLASSES["gebaeude"] + labels[35:50, 35:50] = CLASSES["gebaeude"] + assert len(building_polygons(labels)) == 2 diff --git a/tools/extract/tests/test_heights.py b/tools/extract/tests/test_heights.py new file mode 100644 index 0000000..ea9a7d7 --- /dev/null +++ b/tools/extract/tests/test_heights.py @@ -0,0 +1,54 @@ +from shapely.geometry import Polygon +from leonida.heights import DISTRICT_RANGES, assign_heights + +SQUARE = Polygon([(0, 0), (20, 0), (20, 20), (0, 20)]) +FAR = Polygon([(500, 500), (520, 500), (520, 520), (500, 520)]) + +def test_landmark_height_wins(): + res = assign_heights([SQUARE], landmarks=[{"xyz": [10, 10, 0], "height": 180.0}], + districts=[], seed=1) + assert res[0]["height"] == 180.0 + assert res[0]["source"] == "skyline" + +def test_district_range_is_used_without_landmark(): + res = assign_heights([SQUARE], landmarks=[], + districts=[{"name": "Downtown", "x": 5, "y": 5}], seed=1) + lo, hi = DISTRICT_RANGES["Downtown"] + assert lo <= res[0]["height"] <= hi + assert res[0]["source"] == "district" + +def test_same_seed_gives_same_height(): + args = dict(landmarks=[], districts=[{"name": "Downtown", "x": 5, "y": 5}], seed=42) + assert assign_heights([SQUARE], **args)[0]["height"] == assign_heights([SQUARE], **args)[0]["height"] + +def test_unknown_district_falls_back_to_default(): + res = assign_heights([FAR], landmarks=[], districts=[], seed=3) + lo, hi = DISTRICT_RANGES["default"] + assert lo <= res[0]["height"] <= hi + +from shapely.geometry import LineString, Point +from leonida.heights import PALM_KERB, palm_points + +def test_no_palm_stands_on_the_carriageway(): + """Der Versatz muss vom Fahrbahnrand her gedacht werden, nicht von einem festen Wert. + + Am echten Mosaik hat die Fahrbahn an den Mittellinienpunkten p50 5,4 m, p90 11,7 m + und p95 14,0 m Halbbreite. Ein fester Versatz von 4 m setzt 60 % der Palmen auf die + Fahrbahn, einer von 11 m immer noch 12,5 % – und auf der Hälfte der Straßen stünde + die Palme dann im Nachbargebäude. + """ + road = Polygon([(0, -9), (200, -9), (200, 9), (0, 9)]) # 18 m breite Fahrbahn + line = LineString([(0, 0), (200, 0)]) + sand = Polygon([(0, 20), (200, 20), (200, 40), (0, 40)]) + pts = palm_points([sand], [], [line], [], [road], seed=1) + assert pts + assert all(road.distance(Point(x, y)) >= PALM_KERB for x, y in pts) + +def test_palms_do_not_pile_up_where_centerlines_meet(): + """Mittellinien sind an Kreuzungen getrennt (Task 4). Ein Start bei Position 0 setzt + auf zwei gegenüberliegenden Ästen zwei Palmen in genau denselben Punkt.""" + sand = Polygon([(-80, 20), (80, 20), (80, 40), (-80, 40)]) + arms = [LineString([(0, 0), (60, 0)]), LineString([(0, 0), (-60, 0)])] + pts = palm_points([sand], [], arms, [], [], seed=1) + assert pts + assert len(set(pts)) == len(pts)