Aislo/B04_PreProcess/B04_PreProcess_Engine_Quality.py

"""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 넓이
저장 = `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 = 1
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 _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
    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)