569 lines
25 KiB
Python
569 lines
25 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) -> 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.
|
|
"""
|
|
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)
|
|
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."""
|
|
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)))
|
|
|
|
|
|
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)))))
|