"""B04 서피스 품질 지표 넷 — 기준 없이 재는 값(762-2 · _분석_599 정의 · _분석_601-1 3-3). ① 튐 넓이 % — 5 × 5 칸 중간값 대비 0.5 m 넘는 칸(8 방향 묶음 = 곳 수) ② 1 m 등고선 꺾임 /km — 2 m 마다 다시 찍어 30° 넘게 꺾인 곳 ③ 땅 점 빈 곳 % — 2 m 안 지면점 없는 칸 · 그중 기울기 > 15° 칸 ④ 필터 차 % — 같은 표현 · 짝 필터 모델과 |차| > 0.5 m 넓이 ⑤ 메운 넓이 m² · % — 외곽 안 구멍을 막 펴기로 메운 칸(766-2 · 같은 필터 dtm `filled_mask`) 저장 = `models/quality_{필터}_{표현}.json`(파생 저장본 — 모델 · 짝 모델 지문이 다를 때만 다시 셈). """ import json import logging from pathlib import Path from typing import Any import numpy as np from scipy import ndimage from skimage.measure import find_contours from common_util.common_util_derived_store import file_stamp from common_util.common_util_json import atomic_write_json logger = logging.getLogger(__name__) QUALITY_VERSION = 2 SPIKE_M = 0.5 SPIKE_WINDOW = 5 CONTOUR_INTERVAL_M = 1.0 KINK_STEP_M = 2.0 KINK_DEG = 30.0 GAP_M = 2.0 GAP_STEEP_DEG = 15.0 DIFF_M = 0.5 GRID_METHODS = ("dtm", "implicit") # 필터 차 짝 — 업체 분류 ↔ 점진 TIN(다듬기 판끼리) · 없으면 다듬기 전 ↔ 뒤 _PAIRS = { "classification": ("ptd", "classification_refine"), "ptd": ("classification", "ptd_refine"), "classification_refine": ("ptd_refine", "classification"), "ptd_refine": ("classification_refine", "ptd"), } def _grid_path(models_dir: Path, filter_key: str, method: str) -> Path: suffix = "_grid" if method == "implicit" else "" return models_dir / f"{method}_{filter_key}{suffix}.npz" def _load_grid(path: Path) -> tuple[np.ndarray, np.ndarray, np.ndarray]: with np.load(path) as data: z = np.asarray(data["z"], dtype=np.float64) valid = np.asarray(data["valid_mask"], dtype=bool) & np.isfinite(z) cell = float(np.asarray(data["resolution"]).ravel()[0]) if "resolution" in data else 1.0 return np.where(valid, z, 0.0), valid, np.asarray(cell) def spike_stats(z: np.ndarray, valid: np.ndarray) -> dict[str, Any]: residual = z - ndimage.median_filter(z, size=SPIKE_WINDOW, mode="nearest") inner = ndimage.binary_erosion(valid, iterations=SPIKE_WINDOW // 2) spike = (np.abs(residual) > SPIKE_M) & inner _, count = ndimage.label(spike, structure=np.ones((3, 3), dtype=bool)) area = int(inner.sum()) return { "spike_pct": round(100.0 * int(spike.sum()) / area, 2) if area else None, "spike_count": int(count), } def contour_kinks(z: np.ndarray, valid: np.ndarray, cell: float) -> dict[str, Any]: if not valid.any(): return {"contour_km": 0.0, "kinks_per_km": None, "bubbles": 0} low, high = float(z[valid].min()), float(z[valid].max()) total, kinks, bubbles = 0.0, 0, 0 step = KINK_STEP_M / cell for level in np.arange(np.ceil(low / CONTOUR_INTERVAL_M) * CONTOUR_INTERVAL_M, high, 1.0): for line in find_contours(z, level, mask=valid): seg = np.hypot(*np.diff(line, axis=0).T) length = float(seg.sum()) if length < step * 2: continue total += length if np.allclose(line[0], line[-1]) and length * cell < 30.0: bubbles += 1 along = np.r_[0.0, np.cumsum(seg)] t = np.arange(0.0, length, step) rows, cols = np.interp(t, along, line[:, 0]), np.interp(t, along, line[:, 1]) angle = np.arctan2(np.diff(rows), np.diff(cols)) turn = np.abs((np.diff(angle) + np.pi) % (2 * np.pi) - np.pi) kinks += int((np.degrees(turn) > KINK_DEG).sum()) km = total * cell / 1000.0 return { "contour_km": round(km, 2), "kinks_per_km": round(kinks / km, 1) if km else None, "bubbles": bubbles, } def gap_stats( z: np.ndarray, valid: np.ndarray, occupied: np.ndarray, cell: float ) -> dict[str, Any]: gap = (ndimage.distance_transform_edt(~occupied) * cell > GAP_M) & valid grad_y, grad_x = np.gradient(z, cell) steep = np.degrees(np.arctan(np.hypot(grad_x, grad_y))) > GAP_STEEP_DEG area = int(valid.sum()) return { "gap_pct": round(100.0 * int(gap.sum()) / area, 2) if area else None, "gap_steep_pct": round(100.0 * int((gap & steep).sum()) / area, 2) if area else None, } def ground_occupancy(xyz: Any, indices: np.ndarray, bounds: np.ndarray, shape: tuple) -> np.ndarray: """지면점이 하나라도 든 칸 — dtm 격자와 같은 칸 나눔(`TerrainContext.grid`).""" rows, cols = shape cell_x = (bounds[0, 1] - bounds[0, 0]) / max(cols - 1, 1) cell_y = (bounds[1, 1] - bounds[1, 0]) / max(rows - 1, 1) occupied = np.zeros(shape, dtype=bool) for start in range(0, len(indices), 2_000_000): points = np.asarray(xyz[indices[start : start + 2_000_000]], dtype=np.float64) gx = np.clip(((points[:, 0] - bounds[0, 0]) / cell_x).astype(np.int64), 0, cols - 1) gy = np.clip(((points[:, 1] - bounds[1, 0]) / cell_y).astype(np.int64), 0, rows - 1) occupied[gy, gx] = True return occupied def diff_pct(z: np.ndarray, valid: np.ndarray, other_path: Path) -> float | None: other_z, other_valid, _ = _load_grid(other_path) if other_z.shape != z.shape: return None both = valid & other_valid area = int(both.sum()) if not area: return None return round(100.0 * int(((np.abs(z - other_z) > DIFF_M) & both).sum()) / area, 2) def filled_stats( models_dir: Path, filter_key: str, valid: np.ndarray, cell: float ) -> dict[str, Any]: with np.load(models_dir / f"dtm_{filter_key}.npz") as data: filled = np.asarray(data["filled_mask"], dtype=bool) if "filled_mask" in data else None if filled is None or filled.shape != valid.shape: return {"filled_m2": None, "filled_pct": None} count, area = int((filled & valid).sum()), int(valid.sum()) return { "filled_m2": round(count * cell * cell, 1), "filled_pct": round(100.0 * count / area, 2) if area else None, } def _pair(models_dir: Path, filter_key: str, method: str) -> Path | None: for other in _PAIRS.get(filter_key, ()): path = _grid_path(models_dir, other, method) if path.is_file(): return path return None def ensure_quality( models_dir: Path, filter_key: str, method: str, occupied: np.ndarray | None = None, ) -> dict[str, Any] | None: """지표 넷을 돌려줌 — 저장본 지문이 맞으면 그대로 · ③ 은 지을 때 넘긴 지면점 칸만.""" if method not in GRID_METHODS: return None grid_path = _grid_path(models_dir, filter_key, method) if not grid_path.is_file(): return None pair = _pair(models_dir, filter_key, method) key = { "version": QUALITY_VERSION, "model": file_stamp(grid_path), "pair": [pair.name, file_stamp(pair)] if pair else None, } path = models_dir / f"quality_{filter_key}_{method}.json" try: stored = json.loads(path.read_text(encoding="utf-8")) except (OSError, ValueError): stored = {} if stored.get("key") == key and occupied is None: return stored z, valid, cell = _load_grid(grid_path) cell = float(cell) same_model = (stored.get("key") or {}).get("model") == key["model"] result: dict[str, Any] = {"key": key, "filter": filter_key, "method": method} if same_model and occupied is None: result.update({k: stored.get(k) for k in ("spike_pct", "spike_count", "contour_km")}) result.update({k: stored.get(k) for k in ("kinks_per_km", "bubbles")}) else: result.update(spike_stats(z, valid)) result.update(contour_kinks(z, valid, cell)) if occupied is not None and occupied.shape == z.shape: result.update(gap_stats(z, valid, occupied, cell)) elif same_model: result.update({k: stored.get(k) for k in ("gap_pct", "gap_steep_pct")}) result["diff_against"] = pair.name if pair else None result["diff_pct"] = diff_pct(z, valid, pair) if pair else None result.update(filled_stats(models_dir, filter_key, valid, cell)) atomic_write_json(path, result) return result def ensure_quality_safely(models_dir: Path, filter_key: str, method: str, **kwargs: Any) -> None: """지을 때 부름 — 실패가 모델 짓기를 막지 않음.""" try: ensure_quality(models_dir, filter_key, method, **kwargs) except Exception as exc: # noqa: BLE001 logger.warning("서피스 품질 지표 실패: %s-%s %s", filter_key, method, exc) def quality_after_build( models_dir: Path, filter_key: str, method: str, context: Any, fresh: bool ) -> None: """파이프라인이 표현 하나를 짓거나 다시 쓴 뒤 부름 — ③ 은 지면점이 손에 있을 때만 셈.""" if method not in GRID_METHODS: return occupied = None if fresh or not (models_dir / f"quality_{filter_key}_{method}.json").is_file(): try: resolution = float(context.config.get("dtm_grid_resolution_meters", 1.0)) shape = context.grid(resolution)[2].shape occupied = ground_occupancy(context.xyz, context.indices(), context.bounds, shape) except Exception as exc: # noqa: BLE001 logger.warning("서피스 품질 지면점 칸 실패: %s %s", filter_key, exc) ensure_quality_safely(models_dir, filter_key, method, occupied=occupied)