knsoandClaude Sonnet 5
73a6df7f08
_sides() und _place_rows() nehmen die ersten vier Ringpunkte eines Rechtecks als seine Ecken. GEOS setzt nach dem Puffern je nach Plattform zusätzliche Punkte mitten auf die Kanten (macOS 0, Windows 2), und dann messen beide die falschen Kanten: ergänzte Gebäude bekamen unter entfernten Beschriftungen eine falsche Größe und eine falsche Richtung. simplify(0) entfernt genau diese Punkte und lässt die Form unverändert. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
587 lines
26 KiB
Python
587 lines
26 KiB
Python
"""Wiederaufbau des Klassenbilds unter entfernten Kartenbeschriftungen.
|
|
|
|
Spec: docs/superpowers/specs/2026-09-13-beschriftungen-design.md, Abschnitt „Wiederaufbau“.
|
|
Alle Funktionen arbeiten auf dem Klassenbild aus classify() (1 px = 1 m) und verändern es an Ort
|
|
und Stelle. Pixelkoordinaten sind (Zeile, Spalte), Zeilen wachsen nach Süden. Weltkoordinaten
|
|
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
|
|
from .vectors import _skeleton_graph
|
|
|
|
ROAD_ORDER = ("weg", "spekulativ", "strasse") # Malreihenfolge: später gemalt gewinnt
|
|
AREA_CLASSES = ("wasser", "land", "gruen", "sand", "stadt")
|
|
WINDOW_MARGIN = 60 # m um das Hüllrechteck der Maske
|
|
THIN_DEPTH = 2.0 # m: bis zu diesem Abstand vom Maskenrand gilt eine Maske als dünn
|
|
TRUST = 10.0 # m: Mindestabstand vertrauter Mittellinienpunkte zur Maske
|
|
DIRECTION_LENGTH = 15.0 # m Mittellinie für die Richtung eines Straßenendes
|
|
WIDTH_LENGTH = 5.0 # m Mittellinie für die Halbbreite — weiter hinten liegt oft schon die Kreuzung
|
|
MIN_TRUSTED = 5.0 # m: kürzere vertraute Stücke ergeben kein Straßenende
|
|
REACH = 20.0 # m: so weit darf die Maske vor einem Straßenende liegen
|
|
CROSS_THROUGH = 60.0 # m: größter Abstand zwischen zwei Anschlüssen, über die ein Strahl weiterläuft
|
|
EDGE_MARGIN = 8 # m: so nah am Rand des Ausschnitts ist das Skelett unzuverlässig
|
|
MAX_ANGLE = 20.0 # Grad zwischen Richtung eines Endes und der Verbindungslinie
|
|
WIDTH_RATIO = 1.6 # Halbbreiten eines Paares höchstens um diesen Faktor verschieden
|
|
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
|
|
class RebuildStats:
|
|
connected: int = 0 # Paare von Straßenenden und Enden mit Ziel auf der anderen Seite
|
|
extended: int = 0 # geradeaus weitergeführte Enden, die anschließen
|
|
dead_ends: int = 0 # Enden ohne jeden Anschluss; sie werden nicht gemalt
|
|
closed: int = 0 # geschlossene Gebäude
|
|
added: int = 0 # ergänzte Gebäude
|
|
note: str = ""
|
|
|
|
|
|
@dataclass
|
|
class RebuildResult:
|
|
stats: RebuildStats
|
|
buildings: list[list[tuple[float, float]]] = field(default_factory=list) # Ecken (Zeile, Spalte)
|
|
|
|
|
|
@dataclass(frozen=True)
|
|
class RoadEnd:
|
|
row: float
|
|
col: float
|
|
dr: float # Einheitsvektor in Richtung Maske
|
|
dc: float
|
|
half: float # Halbbreite in m
|
|
|
|
|
|
def fill_nearest(labels: np.ndarray, target: np.ndarray, blocked: np.ndarray) -> None:
|
|
"""Jedes Pixel in target übernimmt die Klasse des nächsten Pixels außerhalb von blocked."""
|
|
if not target.any() or blocked.all():
|
|
return
|
|
index = ndi.distance_transform_edt(blocked, return_distances=False, return_indices=True)
|
|
rows, cols = np.nonzero(target)
|
|
labels[rows, cols] = labels[index[0][rows, cols], index[1][rows, cols]]
|
|
|
|
|
|
def is_thin(mask: np.ndarray) -> bool:
|
|
"""Dünn: kein Maskenpixel liegt weiter als THIN_DEPTH vom Maskenrand entfernt."""
|
|
return float(ndi.distance_transform_edt(np.pad(mask, 1)).max()) <= THIN_DEPTH
|
|
|
|
|
|
def rebuild(labels: np.ndarray, mask: np.ndarray, blocked: np.ndarray | None = None,
|
|
add: bool = True) -> RebuildResult:
|
|
"""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, Flächen,
|
|
Gebäude schließen, Gebäude ergänzen. Mit add=False endet der Aufbau nach dem Schließen;
|
|
ergänzt wird dann über add_buildings, sobald alle Masken aufgebaut sind.
|
|
"""
|
|
blocked = mask if blocked is None else blocked
|
|
stats = RebuildStats()
|
|
win = _window(mask, WINDOW_MARGIN)
|
|
L, M, B = labels[win], mask[win], blocked[win]
|
|
edges = (win[0].start == 0, win[0].stop == labels.shape[0], win[1].start == 0, win[1].stop == labels.shape[1])
|
|
painted = np.zeros(M.shape, dtype=bool)
|
|
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)
|
|
_close_buildings(L, M, B, painted, stats)
|
|
corners = _add_buildings(labels, mask, blocked, stats) if add else []
|
|
return RebuildResult(stats, corners)
|
|
|
|
|
|
def add_buildings(labels: np.ndarray, mask: np.ndarray, blocked: np.ndarray | None = None) -> RebuildResult:
|
|
"""Ergänzt Gebäude unter einer schon aufgebauten Maske; die Statistik trägt nur added und note."""
|
|
blocked = mask if blocked is None else blocked
|
|
stats = RebuildStats()
|
|
corners = _add_buildings(labels, mask, blocked, stats)
|
|
return RebuildResult(stats, corners)
|
|
|
|
|
|
# --- Straßen und Wege -----------------------------------------------------------------------
|
|
|
|
def _plan_roads(L, M, B, edges, stats):
|
|
"""Plant alle Striche, höherrangige Klassen zuerst, damit niedrigere an ihnen anschließen.
|
|
|
|
Gemalt wird danach in ROAD_ORDER, damit an Kreuzungen die Straße gewinnt.
|
|
"""
|
|
dist_mask = ndi.distance_transform_edt(~M)
|
|
any_road = np.isin(L, [CLASSES[name] for name in ROAD_ORDER]) & ~B
|
|
centre = np.zeros(M.shape, dtype=bool)
|
|
strokes = {name: [] for name in ROAD_ORDER}
|
|
for name in reversed(ROAD_ORDER):
|
|
road = CLASSES[name]
|
|
visible = (L == road) & ~B
|
|
if not visible.any():
|
|
continue
|
|
# Die Maske zählt beim Skelettieren als Straße: sonst schrumpft die Mittellinie einer
|
|
# abgeschnittenen Straße um ihre Halbbreite und erreicht die 10-m-Zone nie.
|
|
skeleton = skeletonize(visible | M) & ~M
|
|
graph = _skeleton_graph(skeleton)
|
|
width = ndi.distance_transform_edt(visible)
|
|
ends = [e for e in _road_ends(graph, dist_mask, width) if not _near_edge(e, M.shape, edges)]
|
|
pairs, single = _pair_ends(ends, M | any_road)
|
|
for a, b in pairs:
|
|
points = _hermite(a, b)
|
|
strokes[name].append((road, points, np.linspace(a.half, b.half, len(points))))
|
|
_mark_centre(centre, points)
|
|
stats.connected += 1
|
|
zone = np.argwhere(skeleton & visible & (dist_mask < TRUST))
|
|
for end, length, joins in _cast_rays(L, B, dist_mask, road, single, centre, edges):
|
|
if joins:
|
|
count = max(2, int(math.ceil(length * 2)) + 1)
|
|
steps = np.linspace(0.0, length, count)
|
|
points = list(zip(end.row + end.dr * steps, end.col + end.dc * steps))
|
|
stats.extended += 1
|
|
else:
|
|
points = _to_target(end, zone, M, M | any_road)
|
|
if points is None:
|
|
stats.dead_ends += 1 # Sackgassen werden nicht gemalt
|
|
continue
|
|
stats.connected += 1
|
|
strokes[name].append((road, points, [end.half] * len(points)))
|
|
_mark_centre(centre, points)
|
|
return [stroke for name in ROAD_ORDER for stroke in strokes[name]]
|
|
|
|
|
|
def _to_target(end: RoadEnd, zone, M, passable):
|
|
"""Kurve zur sichtbaren Straße derselben Klasse auf der anderen Seite der Maske, sonst None."""
|
|
best = None
|
|
for r, c in zone:
|
|
vr, vc = r - end.row, c - end.col
|
|
if math.hypot(vr, vc) == 0:
|
|
continue
|
|
angle = _angle((end.dr, end.dc), (vr, vc))
|
|
if angle > MAX_ANGLE or (best is not None and angle >= best[0]):
|
|
continue
|
|
norm = math.hypot(vr, vc)
|
|
target = RoadEnd(float(r), float(c), -vr / norm, -vc / norm, end.half)
|
|
points = _hermite(end, target)
|
|
inside = sum(1 for pr, pc in points if M[int(round(pr)), int(round(pc))])
|
|
if inside >= TARGET_MIN_INSIDE and _corridor(passable, points):
|
|
best = (angle, points)
|
|
return None if best is None else best[1]
|
|
|
|
|
|
def _road_ends(graph, dist_mask, width) -> list[RoadEnd]:
|
|
ends = []
|
|
for node in sorted(graph.nodes):
|
|
if dist_mask[node] < TRUST:
|
|
continue
|
|
inward = [n for n in graph.neighbors(node) if dist_mask[n] < TRUST]
|
|
if not inward:
|
|
continue
|
|
trace = _trace(graph, node, inward, dist_mask)
|
|
if _length(trace) < MIN_TRUSTED:
|
|
continue
|
|
dr, dc = node[0] - trace[-1][0], node[1] - trace[-1][1]
|
|
norm = math.hypot(dr, dc)
|
|
near = [p for k, p in enumerate(trace) if _length(trace[:k + 1]) <= WIDTH_LENGTH]
|
|
half = float(np.median([width[p] for p in near]))
|
|
end = RoadEnd(float(node[0]), float(node[1]), dr / norm, dc / norm, half)
|
|
if not _reaches(dist_mask, end):
|
|
continue
|
|
if not any(math.hypot(e.row - end.row, e.col - end.col) <= 3 and e.dr * end.dr + e.dc * end.dc > 0.9
|
|
for e in ends):
|
|
ends.append(end)
|
|
return ends
|
|
|
|
|
|
def _trace(graph, start, inward, dist_mask):
|
|
"""Mittellinie ab start von der Maske weg, an Verzweigungen möglichst geradeaus."""
|
|
path = [start]
|
|
seen = {start, *inward}
|
|
heading = (start[0] - sum(n[0] for n in inward) / len(inward),
|
|
start[1] - sum(n[1] for n in inward) / len(inward))
|
|
while _length(path) < DIRECTION_LENGTH:
|
|
current = path[-1]
|
|
options = [n for n in graph.neighbors(current) if n not in seen and dist_mask[n] >= TRUST]
|
|
if not options:
|
|
break
|
|
if len(path) >= 3:
|
|
ref = path[max(0, len(path) - 6)]
|
|
heading = (current[0] - ref[0], current[1] - ref[1])
|
|
nxt = max(options, key=lambda n: _cos((n[0] - current[0], n[1] - current[1]), heading))
|
|
path.append(nxt)
|
|
seen.add(nxt)
|
|
return path
|
|
|
|
|
|
def _reaches(dist_mask, end: RoadEnd) -> bool:
|
|
"""Liegt innerhalb von REACH vor dem Ende Maske höchstens eine Halbbreite neben dem Strahl?"""
|
|
h, w = dist_mask.shape
|
|
for s in np.arange(0.5, REACH + 0.25, 0.5):
|
|
r, c = int(round(end.row + end.dr * s)), int(round(end.col + end.dc * s))
|
|
if not (0 <= r < h and 0 <= c < w):
|
|
return False
|
|
if dist_mask[r, c] <= end.half:
|
|
return True
|
|
return False
|
|
|
|
|
|
def _pair_ends(ends, passable):
|
|
"""Paare aufsteigend nach Winkelsumme; die Kurve muss über Maske oder sichtbarer Straße laufen."""
|
|
candidates = []
|
|
for i, a in enumerate(ends):
|
|
for j in range(i + 1, len(ends)):
|
|
b = ends[j]
|
|
vr, vc = b.row - a.row, b.col - a.col
|
|
if math.hypot(vr, vc) == 0:
|
|
continue
|
|
angle_a = _angle((a.dr, a.dc), (vr, vc))
|
|
angle_b = _angle((b.dr, b.dc), (-vr, -vc))
|
|
if angle_a > MAX_ANGLE or angle_b > MAX_ANGLE:
|
|
continue
|
|
if max(a.half, b.half) > WIDTH_RATIO * max(min(a.half, b.half), 0.5):
|
|
continue
|
|
candidates.append((angle_a + angle_b, i, j))
|
|
used, pairs = set(), []
|
|
for _, i, j in sorted(candidates):
|
|
if i in used or j in used or not _corridor(passable, _hermite(ends[i], ends[j])):
|
|
continue
|
|
used.update((i, j))
|
|
pairs.append((ends[i], ends[j]))
|
|
return pairs, [e for k, e in enumerate(ends) if k not in used]
|
|
|
|
|
|
def _corridor(passable, points) -> bool:
|
|
h, w = passable.shape
|
|
hits = [0 <= int(round(r)) < h and 0 <= int(round(c)) < w and passable[int(round(r)), int(round(c))]
|
|
for r, c in points]
|
|
return sum(hits) >= CORRIDOR_SHARE * len(hits)
|
|
|
|
|
|
def _hermite(a: RoadEnd, b: RoadEnd):
|
|
dist = math.hypot(b.row - a.row, b.col - a.col)
|
|
t = np.linspace(0.0, 1.0, max(3, 2 * int(math.ceil(dist)) + 1))
|
|
h00, h10 = 2 * t**3 - 3 * t**2 + 1, t**3 - 2 * t**2 + t
|
|
h01, h11 = -2 * t**3 + 3 * t**2, t**3 - t**2
|
|
rows = h00 * a.row + h10 * a.dr * dist + h01 * b.row - h11 * b.dr * dist
|
|
cols = h00 * a.col + h10 * a.dc * dist + h01 * b.col - h11 * b.dc * dist
|
|
return list(zip(rows, cols))
|
|
|
|
|
|
def _near_edge(end: RoadEnd, shape, edges) -> bool:
|
|
"""Liegt das Ende näher als EDGE_MARGIN an einer Fensterkante, die zugleich Rand des Ausschnitts ist?"""
|
|
top, bottom, left, right = edges
|
|
h, w = shape
|
|
return ((top and end.row < EDGE_MARGIN) or (bottom and end.row > h - 1 - EDGE_MARGIN)
|
|
or (left and end.col < EDGE_MARGIN) or (right and end.col > w - 1 - EDGE_MARGIN))
|
|
|
|
|
|
def _cast_rays(L, B, dist_mask, road, ends, centre, edges):
|
|
"""Weiterführen: (Ende, Länge, schließt an) je Ende ohne Partner.
|
|
|
|
„In der Maske“ heißt: Maske höchstens eine Halbbreite neben der Mittellinie, damit auch eine
|
|
seitlich angeschnittene Straße weiterläuft. Ein Anschluss ist eine gemalte Mittellinie in der
|
|
Maske (Kreuzung) oder sichtbare Straße der eigenen Klasse oder strasse höchstens eine Halbbreite
|
|
neben dem Strahl. Der Strahl läuft von Anschluss zu Anschluss weiter, solange der nächste
|
|
höchstens CROSS_THROUGH entfernt ist, und endet am letzten — im Straßenraster geht eine Straße
|
|
an einer Kreuzung meist weiter. Erreicht er in der Maske den Rand des Ausschnitts, schließt er
|
|
dort an: die Straße läuft über die Karte hinaus. Verlässt er die Maske auf anderem Boden, ist
|
|
Schluss; ohne jeden Anschluss ist das Ende eine Sackgasse.
|
|
"""
|
|
top, bottom, left, right = edges
|
|
h, w = dist_mask.shape
|
|
limit = math.hypot(h, w)
|
|
joinable = ((L == road) | (L == CLASSES["strasse"])) & ~B # Buchstabenkonturen zählen nicht
|
|
rays = []
|
|
for end in ends:
|
|
reach = int(math.ceil(end.half))
|
|
entered, joins, outside = None, [], 0.0
|
|
s = 0.0
|
|
while s <= limit:
|
|
r, c = int(round(end.row + end.dr * s)), int(round(end.col + end.dc * s))
|
|
if not (0 <= r < h and 0 <= c < w):
|
|
if entered is not None and ((r < 0 and top) or (r >= h and bottom) or (c < 0 and left) or (c >= w and right)):
|
|
joins.append(s - 0.5) # Rand des Ausschnitts: die Straße geht weiter
|
|
break
|
|
inside = dist_mask[r, c] <= end.half
|
|
if entered is None:
|
|
entered = s if inside else None
|
|
elif inside:
|
|
outside = 0.0
|
|
if centre[r, c] and dist_mask[r, c] == 0:
|
|
joins.append(s)
|
|
elif joinable[max(0, r - reach):r + reach + 1, max(0, c - reach):c + reach + 1].any():
|
|
joins.append(s)
|
|
outside += 0.5
|
|
if outside > 2 * end.half + 2:
|
|
break # fest an sichtbarer Straße angekommen
|
|
else:
|
|
break # anderer Boden: hier ist Schluss
|
|
if joins and s - joins[-1] > CROSS_THROUGH:
|
|
break
|
|
s += 0.5
|
|
if entered is not None:
|
|
rays.append([end, entered, joins[-1] if joins else s, bool(joins)])
|
|
events = []
|
|
for i in range(len(rays)):
|
|
for j in range(i + 1, len(rays)):
|
|
hit = _crossing(rays[i], rays[j])
|
|
if hit is not None:
|
|
events.append((max(hit), i, j, *hit))
|
|
for _, i, j, ti, tj in sorted(events):
|
|
if ti <= rays[i][2] and tj <= rays[j][2]:
|
|
rays[i][2:] = [ti, True]
|
|
rays[j][2:] = [tj, True]
|
|
return [(end, stop, joins) for end, _, stop, joins in rays]
|
|
|
|
|
|
def _crossing(ray_a, ray_b):
|
|
a, a_in, a_stop, _ = ray_a
|
|
b, b_in, b_stop, _ = ray_b
|
|
denom = a.dr * b.dc - a.dc * b.dr
|
|
if abs(denom) < 1e-9:
|
|
return None
|
|
qr, qc = b.row - a.row, b.col - a.col
|
|
ta = (qr * b.dc - qc * b.dr) / denom
|
|
tb = (qr * a.dc - qc * a.dr) / denom
|
|
if a_in <= ta <= a_stop and b_in <= tb <= b_stop:
|
|
return ta, tb
|
|
return None
|
|
|
|
|
|
def _stamp(L, M, painted, road, points, halves):
|
|
strip = np.zeros(M.shape, dtype=bool)
|
|
for (r, c), half in zip(points, halves):
|
|
rr, cc = disk((r, c), max(float(half), 0.5), shape=M.shape)
|
|
strip[rr, cc] = True
|
|
strip &= M
|
|
L[strip] = road
|
|
painted |= strip
|
|
|
|
|
|
def _mark_centre(centre, points):
|
|
h, w = centre.shape
|
|
for r, c in points:
|
|
r, c = int(round(r)), int(round(c))
|
|
centre[max(0, r - 1):min(h, r + 2), max(0, c - 1):min(w, c + 2)] = True
|
|
|
|
|
|
# --- 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."""
|
|
if not target.any():
|
|
return
|
|
source = np.zeros(L.shape, dtype=bool)
|
|
for name in AREA_CLASSES:
|
|
area = (L == CLASSES[name]) & ~B
|
|
comps, count = ndi.label(area)
|
|
if count:
|
|
sizes = ndi.sum(area, comps, np.arange(1, count + 1))
|
|
source |= np.isin(comps, 1 + np.flatnonzero(sizes >= SOURCE_MIN_AREA))
|
|
if not source.any():
|
|
L[target] = CLASSES["stadt"]
|
|
return
|
|
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):
|
|
rows, cols = np.nonzero(mask)
|
|
h, w = mask.shape
|
|
return (slice(max(0, int(rows.min()) - margin), min(h, int(rows.max()) + 1 + 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.
|
|
|
|
Liegen die Pixel auf einer Linie, ist das Hüllrechteck eine Linie; eckige Enden machen daraus
|
|
wieder ein Rechteck mit vier Ecken statt einer Kapsel.
|
|
"""
|
|
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
|
|
# simplify(0) wirft nur Punkte weg, die genau auf einer Kante liegen, und lässt die Form
|
|
# unverändert. Nötig, weil GEOS je nach Plattform solche Punkte setzt (macOS 0, Windows 2) und
|
|
# _sides() wie _place_rows() die ersten vier Ringpunkte als Ecken nehmen: mit Zusatzpunkten
|
|
# messen sie die falschen Kanten und setzen ergänzte Gebäude in falscher Größe und Richtung.
|
|
return envelope.buffer(0.5, cap_style="square", join_style="mitre").simplify(0)
|
|
|
|
|
|
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)))
|
|
|
|
|
|
def _cos(a, b) -> float:
|
|
na, nb = math.hypot(*a), math.hypot(*b)
|
|
return (a[0] * b[0] + a[1] * b[1]) / (na * nb) if na and nb else -1.0
|
|
|
|
|
|
def _angle(a, b) -> float:
|
|
return math.degrees(math.acos(max(-1.0, min(1.0, _cos(a, b)))))
|