From ce93e579363a0cdacb1d0dff90101e2db42fd067 Mon Sep 17 00:00:00 2001 From: KNSO Date: Sun, 13 Sep 2026 18:55:01 +0200 Subject: [PATCH] =?UTF-8?q?feat(extract):=20angeschnittene=20Geb=C3=A4ude?= =?UTF-8?q?=20schlie=C3=9Fen=20und=20verdeckte=20Blockteile=20bebauen?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 5 --- tools/extract/leonida/inpaint.py | 181 +++++++++++++++++++++++++++- tools/extract/tests/test_inpaint.py | 82 ++++++++++++- 2 files changed, 258 insertions(+), 5 deletions(-) diff --git a/tools/extract/leonida/inpaint.py b/tools/extract/leonida/inpaint.py index 752039d..9119bf0 100644 --- a/tools/extract/leonida/inpaint.py +++ b/tools/extract/leonida/inpaint.py @@ -8,11 +8,14 @@ kennt dieses Modul nicht; die Umrechnung macht annotations.py. from __future__ import annotations import math +import warnings from dataclasses import dataclass, field import numpy as np from scipy import ndimage as ndi +from shapely.geometry import MultiPoint, Polygon from skimage.draw import disk +from skimage.draw import polygon as fill_polygon from skimage.morphology import skeletonize from .classify import CLASSES @@ -34,6 +37,18 @@ WIDTH_RATIO = 1.6 # Halbbreiten eines Paares höchstens um diesen Faktor CORRIDOR_SHARE = 0.9 # Anteil der Kurvenpunkte, die auf Maske oder sichtbarer Straße liegen müssen TARGET_MIN_INSIDE = 3 # Kurvenpunkte in der Maske, ab denen ein Ziel auf der anderen Seite liegt SOURCE_MIN_AREA = 50 # m²: kleinere Flächenstücke sind Kantenglättung und taugen nicht als Füllquelle +TEMPLATE_MIN_AREA = 100 # m²: Gebäude, die als Vorlage und als Nachbarn zählen +NEIGHBOUR_RADIUS = 150 # m um die Maske für Gebäudevorlage und Nachbarschaftsdichte +MIN_NEIGHBOURS = 3 +SHORT_SIDE = (6.0, 40.0) # m, Grenzen der kurzen Seite ergänzter Gebäude +LONG_SIDE_MAX = 60.0 # m +SETBACK = (2.0, 15.0) # m, Grenzen des Straßenabstands +GAP = 3.0 # m zwischen Gebäuden +VISIBLE_SHARE = 0.25 # sichtbarer Blockteil zählt ab 25 % des Blocks … +VISIBLE_AREA = 400 # … und 400 m² +PIECE_MIN_AREA = 12 # m²: kleinere Gebäudestücke werden nicht geschlossen +CONTACT_MIN = 3 # Pixel direkt an der Maske, ab denen ein Gebäudestück angeschnitten ist +EIGHT = np.ones((3, 3), dtype=bool) @dataclass @@ -79,7 +94,8 @@ def rebuild(labels: np.ndarray, mask: np.ndarray, blocked: np.ndarray | None = N """Baut das Klassenbild unter einer zusammenhängenden Maske wieder auf. blocked sind alle Pixel, die nicht als Quelle taugen — alle Beschriftungsmasken; ohne - Angabe nur diese Maske. Reihenfolge wie in der Spec: Straßen und Wege, dann Flächen. + Angabe nur diese Maske. Reihenfolge wie in der Spec: Straßen und Wege, Flächen, + Gebäude schließen, Gebäude ergänzen. """ blocked = mask if blocked is None else blocked stats = RebuildStats() @@ -90,7 +106,9 @@ def rebuild(labels: np.ndarray, mask: np.ndarray, blocked: np.ndarray | None = N for road, points, halves in _plan_roads(L, M, B, edges, stats): _stamp(L, M, painted, road, points, halves) _fill_areas(L, M & ~painted, B) - return RebuildResult(stats) + _close_buildings(L, M, B, painted, stats) + corners = _add_buildings(labels, mask, blocked, stats) + return RebuildResult(stats, corners) # --- Straßen und Wege ----------------------------------------------------------------------- @@ -354,7 +372,7 @@ def _mark_centre(centre, points): centre[max(0, r - 1):min(h, r + 2), max(0, c - 1):min(w, c + 2)] = True -# --- Flächen ------------------------------------------------------------------------------------- +# --- Flächen und Gebäude -------------------------------------------------------------------- def _fill_areas(L, target, B): """Nächste sichtbare Fläche; Stücke unter SOURCE_MIN_AREA (Kantenglättung) zählen nicht als Quelle.""" @@ -373,6 +391,140 @@ def _fill_areas(L, target, B): fill_nearest(L, target, ~source) +def _close_buildings(L, M, B, painted, stats): + building = CLASSES["gebaeude"] + pieces = (L == building) & ~B + comps, count = ndi.label(pieces, structure=EIGHT) + if count == 0: + return + index = np.arange(1, count + 1) + sizes = ndi.sum(pieces, comps, index) + contact = ndi.sum(ndi.binary_dilation(M, structure=EIGHT) & pieces, comps, index) + for k, sl in enumerate(ndi.find_objects(comps), start=1): + if sizes[k - 1] < PIECE_MIN_AREA or contact[k - 1] < CONTACT_MIN: + continue + rows, cols = np.nonzero(comps[sl] == k) + rect = _rectangle(rows + sl[0].start, cols + sl[1].start) + rr, cc = _raster(rect, M.shape) + inside = M[rr, cc] & ~painted[rr, cc] & (L[rr, cc] != building) + if inside.any(): + L[rr[inside], cc[inside]] = building + stats.closed += 1 + + +def _add_buildings(labels, mask, blocked, stats): + """Ergänzt Gebäude im verdeckten Teil der Blöcke; alles im Nachbarschaftsfenster (150 m).""" + building, urban, water = CLASSES["gebaeude"], CLASSES["stadt"], CLASSES["wasser"] + roads = [CLASSES[name] for name in ROAD_ORDER] + hood = _window(mask, NEIGHBOUR_RADIUS + 1) + L, M, B = labels[hood], mask[hood], blocked[hood] + + near = ndi.distance_transform_edt(~M) <= NEIGHBOUR_RADIUS + pieces = (L == building) & ~B + comps, count = ndi.label(pieces, structure=EIGHT) + index = np.arange(1, count + 1) + sizes = ndi.sum(pieces, comps, index) if count else [] + reach = ndi.maximum(near, comps, index) if count else [] + chosen = [k for k in range(1, count + 1) if sizes[k - 1] >= TEMPLATE_MIN_AREA and reach[k - 1]] + if len(chosen) < MIN_NEIGHBOURS: + stats.note = "zu wenig Nachbarn" + return [] + shorts, longs = [], [] + slices = ndi.find_objects(comps) + for k in chosen: + sl = slices[k - 1] + rows, cols = np.nonzero(comps[sl] == k) + short, long_ = _sides(_rectangle(rows + sl[0].start, cols + sl[1].start)) + shorts.append(short) + longs.append(long_) + short = float(np.clip(np.median(shorts), *SHORT_SIDE)) + long_ = float(np.clip(np.median(longs), short, LONG_SIDE_MAX)) + road = np.isin(L, roads) + road_dist = ndi.distance_transform_edt(~road) + neighbours = np.isin(comps, chosen) + outline = neighbours & ~ndi.binary_erosion(neighbours) + setback = float(np.clip(np.median(road_dist[outline]), *SETBACK)) + ground = ~B & near & ~np.isin(L, roads + [water]) + density_hood = float(((L == building) & ~B & near).sum() / max(int(ground.sum()), 1)) + + forbidden = ndi.distance_transform_edt(L != building) < GAP + buildable = (L == urban) & M & (road_dist >= setback) & ~forbidden + blocks, _ = ndi.label(~np.isin(L, roads + [water])) + placed_all = [] + for k, sl in enumerate(ndi.find_objects(blocks), start=1): + block = blocks[sl] == k + hidden = block & M[sl] + if not (hidden & buildable[sl]).any(): + continue + visible = block & ~B[sl] + if visible.sum() >= VISIBLE_SHARE * block.sum() and visible.sum() >= VISIBLE_AREA: + density = ((L[sl] == building) & visible).sum() / visible.sum() + else: + density = density_hood + target = density * hidden.sum() - ((L[sl] == building) & hidden).sum() + if target < 0.5 * short * long_: + continue + rows, cols = np.nonzero(block) + placed = _place_rows(L, buildable, forbidden, rows + sl[0].start, cols + sl[1].start, + road, short, long_, setback, target) + stats.added += len(placed) + placed_all += [[(r + hood[0].start, c + hood[1].start) for r, c in corners] for corners in placed] + return placed_all + + +def _place_rows(L, buildable, forbidden, rows, cols, road, short, long_, setback, target): + """Rechtecke reihenweise von der Straßenseite des Blocks her, bis target erreicht ist.""" + corners = list(_rectangle(rows, cols).exterior.coords)[:4] # (Spalte, Zeile) + p0 = np.array(corners[0]) + e1, e2 = np.array(corners[1]) - p0, np.array(corners[3]) - p0 + if np.hypot(*e1) < np.hypot(*e2): + e1, e2 = e2, e1 + U, V = float(np.hypot(*e1)), float(np.hypot(*e2)) + u, v = e1 / U, e2 / V + origin = p0 + + reach = setback + 2 + rr, cc = np.nonzero(road) + pu = (cc - origin[0]) * u[0] + (rr - origin[1]) * u[1] + pv = (cc - origin[0]) * v[0] + (rr - origin[1]) * v[1] + along = (pu >= 0) & (pu <= U) + near_origin_side = int((along & (pv < 0) & (pv >= -reach)).sum()) + near_far_side = int((along & (pv > V) & (pv <= V + reach)).sum()) + if near_far_side > near_origin_side: + origin, v = origin + e2, -v + + free = buildable[rows, cols] + if not free.any(): + return [] + fu = (cols[free] - origin[0]) * u[0] + (rows[free] - origin[1]) * u[1] + fv = (cols[free] - origin[0]) * v[0] + (rows[free] - origin[1]) * v[1] + u_start, u_end = max(0.0, float(fu.min()) - long_), min(U, float(fu.max())) + v_pos = max(setback, float(fv.min()) - short) + v_end = min(V, float(fv.max())) + + placed, area = [], 0.0 + while v_pos <= v_end and area < target: + u_pos, used = u_start, False + while u_pos <= u_end and area < target: + quad = [origin + u * du + v * dv for du, dv in + ((u_pos, v_pos), (u_pos + long_, v_pos), (u_pos + long_, v_pos + short), (u_pos, v_pos + short))] + polygon = Polygon([tuple(p) for p in quad]) + qr, qc = _raster(polygon, L.shape) + if qr.size and buildable[qr, qc].all() and not forbidden[qr, qc].any(): + L[qr, qc] = CLASSES["gebaeude"] + buildable[qr, qc] = False + gr, gc = _raster(polygon.buffer(GAP, join_style="mitre"), L.shape) + forbidden[gr, gc] = True + placed.append([(float(p[1]), float(p[0])) for p in quad]) + area += short * long_ + used = True + u_pos += long_ + GAP + else: + u_pos += 1.0 + v_pos += short + GAP if used else 1.0 + return placed + + # --- Hilfen ----------------------------------------------------------------------------------- def _window(mask, margin): @@ -382,6 +534,27 @@ def _window(mask, margin): slice(max(0, int(cols.min()) - margin), min(w, int(cols.max()) + 1 + margin))) +def _rectangle(rows, cols) -> Polygon: + """Gedrehtes Hüllrechteck um Pixel (Zeilen, Spalten) als Polygon in (Spalte, Zeile), mit Pixelkante.""" + points = MultiPoint(np.column_stack([cols, rows]).astype(float)) + with warnings.catch_warnings(): + # Siehe buildings.regularize(): harmlose GEOS-Warnung bei achsenparallelen Formen. + warnings.filterwarnings("ignore", message=r".*oriented_envelope.*", category=RuntimeWarning) + envelope = points.minimum_rotated_rectangle + return envelope.buffer(0.5, join_style="mitre") + + +def _sides(rect: Polygon) -> tuple[float, float]: + corners = list(rect.exterior.coords)[:4] + edges = sorted(math.dist(corners[i - 1], corners[i]) for i in range(4)) + return edges[0], edges[-1] + + +def _raster(polygon: Polygon, shape): + xs, ys = polygon.exterior.coords.xy + return fill_polygon(np.asarray(ys), np.asarray(xs), shape=shape) + + def _length(path) -> float: return sum(math.dist(path[i - 1], path[i]) for i in range(1, len(path))) @@ -392,4 +565,4 @@ def _cos(a, b) -> float: def _angle(a, b) -> float: - return math.degrees(math.acos(max(-1.0, min(1.0, _cos(a, b))))) \ No newline at end of file + return math.degrees(math.acos(max(-1.0, min(1.0, _cos(a, b))))) diff --git a/tools/extract/tests/test_inpaint.py b/tools/extract/tests/test_inpaint.py index c3e9f11..804bd7e 100644 --- a/tools/extract/tests/test_inpaint.py +++ b/tools/extract/tests/test_inpaint.py @@ -161,4 +161,84 @@ def test_a_road_running_off_the_map_under_text_is_kept_to_the_edge(): scribble(labels, mask) result = rebuild(labels, mask) assert result.stats.extended == 1 - assert (labels[0:30, 142:158] == ROAD).all() \ No newline at end of file + assert (labels[0:30, 142:158] == ROAD).all() + + +def test_a_bitten_building_is_closed_but_not_extended(): + labels = city() + labels[50:90, 50:110] = BUILDING + mask = box_mask(labels.shape, slice(70, 110), slice(90, 150)) + scribble(labels, mask) + result = rebuild(labels, mask) + assert result.stats.closed == 1 + assert (labels[70:90, 90:110] == BUILDING).all() + assert not (labels[92:110, 90:150] == BUILDING).any() + assert not (labels[70:90, 112:150] == BUILDING).any() + + +def block_city(): + """Ein Block zwischen vier Straßen, links sichtbar bebaut, rechts verdeckt.""" + labels = city(300, 400) + labels[20:40, :] = ROAD + labels[260:280, :] = ROAD + labels[:, 20:40] = ROAD + labels[:, 360:380] = ROAD + for top in (60, 100, 140, 180): + for left in (60, 120): + labels[top:top + 20, left:left + 30] = BUILDING + mask = box_mask(labels.shape, slice(60, 241), slice(220, 341)) + return labels, mask + + +def test_a_hidden_block_part_gets_the_density_of_the_visible_part(): + labels, mask = block_city() + visible_density = 8 * 600 / (220 * 320 - mask.sum()) + scribble(labels, mask) + result = rebuild(labels, mask) + added = (labels[mask] == BUILDING).sum() + assert result.stats.added == len(result.buildings) >= 3 + assert abs(added - visible_density * mask.sum()) <= 600 + 2 * (20 + 30) # ein Vorlage-Rechteck genau + + +def test_added_buildings_keep_setback_gaps_and_stay_in_the_hidden_part(): + labels, mask = block_city() + scribble(labels, mask) + result = rebuild(labels, mask) + new = (labels == BUILDING) & mask + assert not (new & ~mask).any() + road_distance = ndi.distance_transform_edt(~np.isin(labels, [ROAD, PATH])) + assert road_distance[new].min() >= 15 # Straßenabstand der Nachbarn, auf 15 m begrenzt + pieces, count = ndi.label(new) + assert count == len(result.buildings) + for k in range(1, count + 1): + others = new & (pieces != k) + assert ndi.distance_transform_edt(~others)[pieces == k].min() >= 3 + + +def test_no_buildings_are_added_in_a_park(): + labels, mask = block_city() + labels[40:260, 200:360] = PARK + scribble(labels, mask) + result = rebuild(labels, mask) + assert result.stats.added == 0 + + +def test_without_three_neighbours_nothing_is_added(): + labels, mask = block_city() + labels[labels == BUILDING] = STADT + labels[60:80, 60:90] = BUILDING + labels[100:120, 60:90] = BUILDING + scribble(labels, mask) + result = rebuild(labels, mask) + assert result.stats.added == 0 + assert result.stats.note == "zu wenig Nachbarn" + + +def test_rebuild_is_deterministic(): + first, mask = block_city() + second = first.copy() + scribble(first, mask) + scribble(second, mask) + a, b = rebuild(first, mask), rebuild(second, mask) + assert (first == second).all() + assert a.buildings == b.buildings