101 lines
3.8 KiB
Python
101 lines
3.8 KiB
Python
import numpy as np
|
|
from scipy import ndimage as ndi
|
|
from shapely.geometry import LineString, Polygon
|
|
from skimage.measure import find_contours
|
|
from skimage.morphology import closing, disk, skeletonize
|
|
|
|
from .classify import CLASSES
|
|
from .coords import SLICE
|
|
|
|
ROAD_CLASSES = ("strasse", "weg", "spekulativ")
|
|
|
|
|
|
def _to_world(contour: np.ndarray) -> list[tuple[float, float]]:
|
|
"""Konturpunkte (Zeile, Spalte) → Weltkoordinaten in Metern."""
|
|
return [(SLICE["x0"] + c, SLICE["y1"] - r) for r, c in contour]
|
|
|
|
|
|
def mask_polygons(mask: np.ndarray, min_area: float, simplify: float) -> list[Polygon]:
|
|
padded = np.pad(mask, 1)
|
|
outers, holes = [], []
|
|
for contour in find_contours(padded.astype(float), 0.5):
|
|
local = contour - 1
|
|
ring = _to_world(local)
|
|
if len(ring) < 4:
|
|
continue
|
|
poly = Polygon(ring).buffer(0)
|
|
if poly.is_empty or poly.area < min_area:
|
|
continue
|
|
(outers if _signed_area(local) > 0 else holes).append(poly)
|
|
result = []
|
|
for outer in outers:
|
|
inner = [h.exterior.coords for h in holes if outer.contains(h)]
|
|
result.append(Polygon(outer.exterior.coords, inner).buffer(0).simplify(simplify))
|
|
return [p for p in result if p.geom_type == "Polygon" and p.area >= min_area]
|
|
|
|
|
|
def _signed_area(contour: np.ndarray) -> float:
|
|
"""Shoelace-Summe in Pixelkoordinaten (Zeile, Spalte).
|
|
|
|
find_contours liefert Außenkonturen und Löcher mit entgegengesetztem Umlaufsinn;
|
|
das Vorzeichen trennt sie zuverlässig. Eine Stichprobe der Maske am Mittelpunkt der
|
|
Kontur (naheliegende Alternative) versagt bei konzentrischen Formen wie einer Insel
|
|
mit Lagune: der Mittelpunkt der äußeren Kontur liegt dann selbst in der Lagune.
|
|
"""
|
|
row, col = contour[:, 0], contour[:, 1]
|
|
return float(0.5 * np.sum(col * np.roll(row, -1) - np.roll(col, -1) * row))
|
|
|
|
|
|
def land_polygons(labels: np.ndarray) -> list[Polygon]:
|
|
land = labels != CLASSES["wasser"]
|
|
land = ndi.binary_closing(ndi.binary_opening(land, disk(2)), disk(2))
|
|
return mask_polygons(land, min_area=50, simplify=1.0)
|
|
|
|
|
|
def road_mask(labels: np.ndarray) -> np.ndarray:
|
|
mask = np.zeros(labels.shape, dtype=bool)
|
|
for name in ROAD_CLASSES:
|
|
mask |= labels == CLASSES[name]
|
|
# Schließen überbrückt die eingezeichneten Straßennamen.
|
|
return closing(mask, disk(3))
|
|
|
|
|
|
# Rand in Pixeln (= Metern) um die Straßenfläche zur Brückenprüfung: die Straßenfläche
|
|
# selbst ist per Definition nie "wasser", daher muss das Fenster über ihren Rand
|
|
# hinausschauen, um das Wasser zu sehen, das sie überquert.
|
|
BRIDGE_MARGIN_PX = 3
|
|
|
|
|
|
def road_polygons(labels: np.ndarray) -> list[tuple[Polygon, bool]]:
|
|
water = labels == CLASSES["wasser"]
|
|
polys = mask_polygons(road_mask(labels), min_area=20, simplify=0.8)
|
|
out = []
|
|
for poly in polys:
|
|
r0, c0, r1, c1 = _bbox_px(poly)
|
|
r0 = max(r0 - BRIDGE_MARGIN_PX, 0)
|
|
c0 = max(c0 - BRIDGE_MARGIN_PX, 0)
|
|
r1 = min(r1 + BRIDGE_MARGIN_PX, water.shape[0] - 1)
|
|
c1 = min(c1 + BRIDGE_MARGIN_PX, water.shape[1] - 1)
|
|
window = water[r0:r1 + 1, c0:c1 + 1]
|
|
out.append((poly, bool(window.mean() > 0.3)))
|
|
return out
|
|
|
|
|
|
def _bbox_px(poly: Polygon) -> tuple[int, int, int, int]:
|
|
minx, miny, maxx, maxy = poly.bounds
|
|
return (int(SLICE["y1"] - maxy), int(minx - SLICE["x0"]),
|
|
int(SLICE["y1"] - miny), int(maxx - SLICE["x0"]))
|
|
|
|
|
|
def centerlines(labels: np.ndarray) -> list[LineString]:
|
|
"""Skelett der Straßenflächen, nur zum Platzieren von Laternen und Palmen."""
|
|
skeleton = skeletonize(road_mask(labels))
|
|
lines = []
|
|
for contour in find_contours(skeleton.astype(float), 0.5):
|
|
pts = _to_world(contour)
|
|
if len(pts) >= 2:
|
|
line = LineString(pts).simplify(1.5)
|
|
if line.length >= 15:
|
|
lines.append(line)
|
|
return lines
|