feat(extract): angeschnittene Gebäude schließen und verdeckte Blockteile bebauen

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

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)))))
+    return math.degrees(math.acos(max(-1.0, min(1.0, _cos(a, 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()
+    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