"""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)))))