"""B05 지형 스켈레톤(주/지 능선·계곡) 추출.

수문학적 정의: D8 흐름누적이 임계값 이상인 셀=계곡, DEM 반전 시 능선.
누적값 크기의 2차 임계값으로 주/지를 나눈다. whitebox 우선, 실패 시 numpy
D8 폴백. 산출 polyline은 비용면과 동일 좌표계이며 B05_wf2_Route/route에 캐시된다.
"""

import json
import math
import tempfile
import time
from pathlib import Path
from typing import Any

import numpy as np

from B05_wf2_Route.B05_wf2_Route_Engine_Solver import (
    _MODELS_SUBDIR,
    _ROUTE_CACHE_SUBDIR,
    _cost_surface_signature,
    _load_or_build_cost_surface,
)
from config.config_system import (
    SKELETON_MAIN_RIDGE_ACC_THRESHOLD_CELLS,
    SKELETON_MAIN_VALLEY_ACC_THRESHOLD_CELLS,
    SKELETON_RIDGE_ACC_THRESHOLD_CELLS,
    SKELETON_VALLEY_ACC_THRESHOLD_CELLS,
)

SKELETON_CLASSES = ("main_ridge", "minor_ridge", "main_valley", "minor_valley")

_D8_OFFSETS = [
    (-1, -1),
    (-1, 0),
    (-1, 1),
    (0, -1),
    (0, 1),
    (1, -1),
    (1, 0),
    (1, 1),
]


def _default_thresholds() -> dict[str, float]:
    return {
        "valley_acc": float(SKELETON_VALLEY_ACC_THRESHOLD_CELLS),
        "main_valley_acc": float(SKELETON_MAIN_VALLEY_ACC_THRESHOLD_CELLS),
        "ridge_acc": float(SKELETON_RIDGE_ACC_THRESHOLD_CELLS),
        "main_ridge_acc": float(SKELETON_MAIN_RIDGE_ACC_THRESHOLD_CELLS),
    }


def d8_flow_accumulation_numpy(z_grid: np.ndarray, valid_mask: np.ndarray) -> np.ndarray:
    """numpy 기반 D8 흐름누적 (whitebox 폴백)."""
    rows, cols = z_grid.shape
    n = rows * cols
    z_flat = z_grid.ravel()
    valid_flat = valid_mask.ravel()

    receiver = np.full(n, -1, dtype=np.int64)
    for idx in range(n):
        if not valid_flat[idx]:
            continue
        r, c = divmod(idx, cols)
        zc = z_flat[idx]
        best_slope = 0.0
        best = -1
        for dr, dc in _D8_OFFSETS:
            nr, nc = r + dr, c + dc
            if not (0 <= nr < rows and 0 <= nc < cols):
                continue
            nidx = nr * cols + nc
            if not valid_flat[nidx]:
                continue
            drop = zc - z_flat[nidx]
            if drop <= 0:
                continue
            dist = math.sqrt(2.0) if (dr != 0 and dc != 0) else 1.0
            slope = drop / dist
            if slope > best_slope:
                best_slope = slope
                best = nidx
        receiver[idx] = best

    acc = np.where(valid_flat, 1.0, 0.0)
    order = np.argsort(-z_flat, kind="stable")
    for idx in order:
        recv = receiver[idx]
        if recv >= 0 and valid_flat[idx]:
            acc[recv] += acc[idx]
    return acc.reshape(rows, cols)


def _d8_flow_accumulation_whitebox(
    z_grid: np.ndarray,
    x_coords: np.ndarray,
    y_coords: np.ndarray,
    valid_mask: np.ndarray,
    grid_res: float,
) -> np.ndarray | None:
    """WhiteboxTools로 D8 흐름누적을 계산한다 (실패 시 None → numpy 폴백)."""
    try:
        import rasterio
        from rasterio.transform import from_origin
        from whitebox import WhiteboxTools
    except Exception:
        return None

    nodata = -9999.0
    rows, cols = z_grid.shape
    z_out = np.where(valid_mask, z_grid, nodata).astype(np.float32)[::-1, :]
    transform = from_origin(
        float(x_coords[0]) - grid_res / 2.0,
        float(y_coords[-1]) + grid_res / 2.0,
        grid_res,
        grid_res,
    )
    try:
        with tempfile.TemporaryDirectory(prefix="wbt_skel_") as tmp:
            tmp_path = Path(tmp)
            dem_tif = tmp_path / "dem.tif"
            acc_tif = tmp_path / "acc.tif"
            with rasterio.open(
                dem_tif,
                "w",
                driver="GTiff",
                height=rows,
                width=cols,
                count=1,
                dtype="float32",
                nodata=nodata,
                transform=transform,
            ) as dst:
                dst.write(z_out, 1)

            wbt = WhiteboxTools()
            wbt.set_verbose_mode(False)
            wbt.set_working_dir(str(tmp_path))
            if wbt.fill_depressions("dem.tif", "filled.tif") != 0:
                raise RuntimeError("fill_depressions 실패")
            if wbt.d8_flow_accumulation("filled.tif", "acc.tif", out_type="cells") != 0:
                raise RuntimeError("d8_flow_accumulation 실패")

            with rasterio.open(acc_tif) as src:
                acc = src.read(1).astype(np.float64)[::-1, :]
            return np.where(np.isfinite(acc) & (acc > 0) & valid_mask, acc, 0.0)
    except Exception:
        return None


def _flow_accumulation(
    z_grid: np.ndarray,
    x_coords: np.ndarray,
    y_coords: np.ndarray,
    valid_mask: np.ndarray,
    grid_res: float,
) -> np.ndarray:
    acc = _d8_flow_accumulation_whitebox(z_grid, x_coords, y_coords, valid_mask, grid_res)
    if acc is None:
        acc = d8_flow_accumulation_numpy(z_grid, valid_mask)
    return acc


def _trace_polylines(mask: np.ndarray) -> list[list[tuple[int, int]]]:
    """1픽셀 폭 스켈레톤 마스크를 (r, c) polyline 목록으로 변환한다."""
    pixels = set(zip(*np.nonzero(mask)))
    if not pixels:
        return []

    def neighbors(p):
        r, c = p
        return [(r + dr, c + dc) for dr, dc in _D8_OFFSETS if (r + dr, c + dc) in pixels]

    degree = {p: len(neighbors(p)) for p in pixels}
    seeds = [p for p in pixels if degree[p] != 2]
    visited_edges: set = set()
    polylines: list[list[tuple[int, int]]] = []

    def edge_key(a, b):
        return (a, b) if a <= b else (b, a)

    def walk(start, nxt):
        path = [start, nxt]
        visited_edges.add(edge_key(start, nxt))
        prev, curr = start, nxt
        while degree[curr] == 2:
            candidates = [q for q in neighbors(curr) if q != prev]
            if not candidates:
                break
            q = candidates[0]
            if edge_key(curr, q) in visited_edges:
                break
            visited_edges.add(edge_key(curr, q))
            path.append(q)
            prev, curr = curr, q
        return path

    for seed in seeds:
        for nb in neighbors(seed):
            if edge_key(seed, nb) not in visited_edges:
                polylines.append(walk(seed, nb))

    for p in pixels:
        if degree[p] == 2:
            for nb in neighbors(p):
                if edge_key(p, nb) not in visited_edges:
                    polylines.append(walk(p, nb))

    return [pl for pl in polylines if len(pl) >= 2]


def _mask_to_polylines(
    mask: np.ndarray, x_coords: np.ndarray, y_coords: np.ndarray, z_grid: np.ndarray
) -> list[dict[str, Any]]:
    """셀 마스크를 세선화한 뒤 모델좌표 polyline 목록으로 변환한다."""
    if not mask.any():
        return []
    try:
        from skimage.morphology import skeletonize

        skel = skeletonize(mask)
    except Exception:
        skel = mask
    out = []
    for pixel_path in _trace_polylines(skel):
        poly = [
            [float(x_coords[c]), float(y_coords[r]), float(z_grid[r, c])] for r, c in pixel_path
        ]
        out.append({"polyline": poly})
    return out


def extract_skeleton_from_grid(
    x_coords: np.ndarray,
    y_coords: np.ndarray,
    z_grid: np.ndarray,
    valid_mask: np.ndarray,
    grid_res: float,
    thresholds: dict[str, float] | None = None,
    use_whitebox: bool = True,
) -> dict[str, list[dict[str, Any]]]:
    """격자에서 주/지 능선·계곡 polyline을 추출한다 (캐시 없음, 테스트용 공개 API)."""
    th = thresholds or _default_thresholds()

    if use_whitebox:
        acc_valley = _flow_accumulation(z_grid, x_coords, y_coords, valid_mask, grid_res)
        acc_ridge = _flow_accumulation(-z_grid, x_coords, y_coords, valid_mask, grid_res)
    else:
        acc_valley = d8_flow_accumulation_numpy(z_grid, valid_mask)
        acc_ridge = d8_flow_accumulation_numpy(-z_grid, valid_mask)

    valley_mask = valid_mask & (acc_valley >= th["valley_acc"])
    main_valley_mask = valley_mask & (acc_valley >= th["main_valley_acc"])
    minor_valley_mask = valley_mask & ~main_valley_mask

    ridge_mask = valid_mask & ~valley_mask & (acc_ridge >= th["ridge_acc"])
    main_ridge_mask = ridge_mask & (acc_ridge >= th["main_ridge_acc"])
    minor_ridge_mask = ridge_mask & ~main_ridge_mask

    return {
        "main_ridge": _mask_to_polylines(main_ridge_mask, x_coords, y_coords, z_grid),
        "minor_ridge": _mask_to_polylines(minor_ridge_mask, x_coords, y_coords, z_grid),
        "main_valley": _mask_to_polylines(main_valley_mask, x_coords, y_coords, z_grid),
        "minor_valley": _mask_to_polylines(minor_valley_mask, x_coords, y_coords, z_grid),
    }


def _skeleton_signature(models_dir: Path, filter_key: str, method: str, smooth: bool) -> str:
    """소스 모델 + 격자 해상도 + 분류 임계값 서명."""
    th = _default_thresholds()
    th_part = "|".join(f"{k}={v}" for k, v in sorted(th.items()))
    return _cost_surface_signature(models_dir, filter_key, method, smooth) + "|" + th_part


def load_or_build_skeleton(
    project_root: Path, filter_key: str, method: str, smooth: bool
) -> dict[str, Any]:
    """스켈레톤을 캐시에서 로드하거나 새로 계산해 저장한다."""
    project_root = Path(project_root)
    models_dir = project_root / _MODELS_SUBDIR
    cache_dir = project_root / _ROUTE_CACHE_SUBDIR
    suffix = "_smooth" if smooth else ""
    cache_path = cache_dir / f"terrain_skeleton_{filter_key}_{method}{suffix}.json"
    signature = _skeleton_signature(models_dir, filter_key, method, smooth)

    if cache_path.exists():
        try:
            with open(cache_path, encoding="utf-8") as f:
                cached = json.load(f)
            if cached.get("signature") == signature:
                return cached
        except Exception:
            pass

    _t0 = time.time()
    (x_coords, y_coords, z_grid, valid_mask, _dz_dx, _dz_dy, grid_res) = (
        _load_or_build_cost_surface(project_root, models_dir, filter_key, method, smooth)
    )
    z_grid = np.asarray(z_grid, dtype=np.float64)
    valid_mask = np.asarray(valid_mask, dtype=bool)

    skeleton = extract_skeleton_from_grid(
        np.asarray(x_coords), np.asarray(y_coords), z_grid, valid_mask, float(grid_res)
    )
    result: dict[str, Any] = dict(skeleton)
    result["grid_res"] = float(grid_res)
    result["signature"] = signature

    try:
        cache_dir.mkdir(parents=True, exist_ok=True)
        with open(cache_path, "w", encoding="utf-8") as f:
            json.dump(result, f, ensure_ascii=False)
    except Exception:
        pass
    return result