"""실제 지형 DEM으로 DXF TIN 외곽을 확장하는 유틸리티.

동기:
    사용자가 제공한 DXF 평면도의 등고선 범위를 벗어나면 3D 장면에서 지형이
    절벽처럼 끊긴다. 조감도에서 이 구역을 AI 프롬프트로 메우면 사실 기반이
    아니라 AI 상상이라 결과가 가변적이다. 이 모듈은 DXF 범위 바깥에 대해
    실제 공개 DEM(AWS Open Terrain Tiles, terrarium PNG 포맷)을 받아 외곽 링
    메시를 만들고, 기존 TIN과 Z-연속성을 맞춰 이어 붙인다.

기본 소스:
    AWS Open Terrain Tiles (terrarium): API 키 없음, 글로벌, ~30m 수직 정확도.
    URL: https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png
    디코딩: elev_m = (R*256 + G + B/256) - 32768

옵션 소스:
    - 로컬 GeoTIFF: `cache/dem/local.tif` 가 존재하고 rasterio가 설치되어
      있으면 우선 사용. 한국이면 NGII 5m DEM을 여기 놓아두면 정확도 상승.

좌표계:
    - TIN/메시 프레임: 투영 CRS(예: EPSG:5187) - origin (zero-basing)
    - DEM fetch: WGS84 (EPSG:4326)
    - 본 모듈은 투영 CRS와 origin만 받으면 내부에서 pyproj로 변환 처리.

Seam 처리:
    1. 수직 datum 자동 보정: inner 경계 근처 TIN Z와 DEM Z의 평균 차이로
       global offset을 잡아 DEM 전체 Z에서 차감.
    2. Feathering: inner 경계로부터 feather_m 이내 도넛 점은 가장 가까운
       TIN 정점의 Z와 보정된 DEM Z를 선형 블렌드.
"""
from __future__ import annotations

import hashlib
import io
import math
from dataclasses import dataclass
from pathlib import Path
from collections.abc import Callable

import numpy as np
import pyproj
import pyvista as pv
import requests
from PIL import Image
from scipy.spatial import Delaunay, cKDTree


# --- AWS Open Terrain Tiles (terrarium) ---------------------------------
_AWS_TERRARIUM_URL = "https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png"
_HEADERS = {
    "User-Agent": "Mozilla/5.0 (S-CANVAS DEM extender)",
    "Accept": "image/png,image/*;q=0.8",
}


@dataclass
class DemExtendResult:
    mesh: pv.PolyData
    n_points: int
    n_faces: int
    buffer_m: float
    source: str
    vertical_offset_m: float = 0.0
    zoom: int = 0
    grid_step_m: float = 0.0
    feather_m: float = 0.0
    info: str = ""


def _latlon_to_tile(lat: float, lon: float, zoom: int) -> tuple[int, int]:
    lat_rad = math.radians(lat)
    n = 2 ** zoom
    x = int((lon + 180.0) / 360.0 * n)
    y = int((1.0 - math.log(math.tan(lat_rad) + 1.0 / math.cos(lat_rad)) / math.pi) / 2.0 * n)
    return x, y


def _tile_to_latlon(x: int, y: int, zoom: int) -> tuple[float, float]:
    n = 2 ** zoom
    lon = x / n * 360.0 - 180.0
    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * y / n)))
    return math.degrees(lat_rad), lon


def _terrarium_decode(rgb: np.ndarray) -> np.ndarray:
    """(H,W,3) uint8 RGB → (H,W) float elevation in meters."""
    r = rgb[..., 0].astype(np.float64)
    g = rgb[..., 1].astype(np.float64)
    b = rgb[..., 2].astype(np.float64)
    return (r * 256.0 + g + b / 256.0) - 32768.0


def _bbox_cache_key(min_lat: float, min_lon: float, max_lat: float, max_lon: float, zoom: int) -> str:
    raw = f"{min_lat:.6f}_{min_lon:.6f}_{max_lat:.6f}_{max_lon:.6f}_z{zoom}"
    return hashlib.sha1(raw.encode("utf-8")).hexdigest()[:16]


def fetch_terrarium_grid(min_lat: float, min_lon: float, max_lat: float, max_lon: float,
                         zoom: int = 13,
                         cache_dir: str | Path = "cache/dem",
                         log_fn: Callable[[str], None] = print,
                         timeout: float = 15.0) -> tuple[np.ndarray, tuple[float, float, float, float]]:
    """AWS terrarium 타일을 BBOX 범위로 받아 합쳐 (elev_grid, bounds) 반환.

    Returns:
        elev_grid: (H,W) float64, 위→아래로 lat 감소, 좌→우로 lon 증가.
        bounds: (grid_lat_max, grid_lon_min, grid_lat_min, grid_lon_max) — 타일
                경계에 정렬된 실제 격자 범위.
    """
    cache_dir = Path(cache_dir)
    cache_dir.mkdir(parents=True, exist_ok=True)
    cache_key = _bbox_cache_key(min_lat, min_lon, max_lat, max_lon, zoom)
    cache_path = cache_dir / f"terrarium_{cache_key}.png"
    bounds_path = cache_dir / f"terrarium_{cache_key}.bounds.txt"

    x_min, y_min = _latlon_to_tile(max_lat, min_lon, zoom)
    x_max, y_max = _latlon_to_tile(min_lat, max_lon, zoom)
    cols = x_max - x_min + 1
    rows = y_max - y_min + 1

    grid_lat_max, grid_lon_min = _tile_to_latlon(x_min, y_min, zoom)
    grid_lat_min, grid_lon_max = _tile_to_latlon(x_max + 1, y_max + 1, zoom)
    bounds = (grid_lat_max, grid_lon_min, grid_lat_min, grid_lon_max)

    if cache_path.exists() and bounds_path.exists():
        try:
            img = Image.open(cache_path).convert("RGB")
            arr = np.asarray(img, dtype=np.uint8)
            elev = _terrarium_decode(arr)
            log_fn(f"  [DEM] 캐시 사용: {cache_path.name} ({arr.shape[1]}x{arr.shape[0]}px, {cols}x{rows} tiles z{zoom})")
            return elev, bounds
        except Exception as e:
            log_fn(f"  [DEM] 캐시 로드 실패 ({e}), 재다운로드합니다.")

    tile_size = 256
    merged = Image.new("RGB", (cols * tile_size, rows * tile_size), (128, 0, 0))
    ok = 0
    fail = 0
    for ty in range(y_min, y_max + 1):
        for tx in range(x_min, x_max + 1):
            url = _AWS_TERRARIUM_URL.format(z=zoom, x=tx, y=ty)
            try:
                resp = requests.get(url, headers=_HEADERS, timeout=timeout)
                if resp.status_code == 200 and len(resp.content) > 200:
                    tile_img = Image.open(io.BytesIO(resp.content)).convert("RGB")
                    merged.paste(tile_img, ((tx - x_min) * tile_size, (ty - y_min) * tile_size))
                    ok += 1
                else:
                    fail += 1
            except requests.exceptions.RequestException:
                fail += 1

    log_fn(f"  [DEM] 타일: 성공 {ok}장, 실패 {fail}장 (z{zoom}, {cols}x{rows})")
    if ok == 0:
        raise RuntimeError("AWS terrarium 타일을 한 장도 받지 못했습니다. 네트워크를 확인하세요.")

    merged.save(cache_path)
    bounds_path.write_text(
        f"{grid_lat_max} {grid_lon_min} {grid_lat_min} {grid_lon_max} z{zoom}\n",
        encoding="utf-8",
    )

    arr = np.asarray(merged, dtype=np.uint8)
    elev = _terrarium_decode(arr)
    return elev, bounds


def _sample_grid_bilinear(elev: np.ndarray,
                          bounds: tuple[float, float, float, float],
                          lats: np.ndarray, lons: np.ndarray) -> np.ndarray:
    """(H,W) 격자에서 (lats, lons) 점들의 Z를 bilinear 샘플링."""
    grid_lat_max, grid_lon_min, grid_lat_min, grid_lon_max = bounds
    h, w = elev.shape
    fx = (lons - grid_lon_min) / (grid_lon_max - grid_lon_min) * (w - 1)
    fy = (grid_lat_max - lats) / (grid_lat_max - grid_lat_min) * (h - 1)
    fx = np.clip(fx, 0.0, w - 1 - 1e-9)
    fy = np.clip(fy, 0.0, h - 1 - 1e-9)
    x0 = np.floor(fx).astype(np.int64); x1 = x0 + 1
    y0 = np.floor(fy).astype(np.int64); y1 = y0 + 1
    x1 = np.clip(x1, 0, w - 1)
    y1 = np.clip(y1, 0, h - 1)
    tx = fx - x0
    ty = fy - y0
    v00 = elev[y0, x0]; v01 = elev[y0, x1]
    v10 = elev[y1, x0]; v11 = elev[y1, x1]
    v0 = v00 * (1 - tx) + v01 * tx
    v1 = v10 * (1 - tx) + v11 * tx
    return v0 * (1 - ty) + v1 * ty


def _try_local_geotiff(cache_dir: Path) -> tuple[np.ndarray, tuple[float, float, float, float], str] | None:
    """cache/dem/local.tif 가 있으면 읽어 (elev, (lat_max,lon_min,lat_min,lon_max), label) 반환.

    rasterio가 없거나 파일 없으면 None.
    """
    tif = cache_dir / "local.tif"
    if not tif.exists():
        return None
    try:
        import rasterio  # type: ignore
        from rasterio.warp import transform_bounds  # type: ignore
    except ImportError:
        return None
    try:
        with rasterio.open(tif) as ds:
            arr = ds.read(1).astype(np.float64)
            # nodata 처리
            nodata = ds.nodata
            if nodata is not None:
                arr = np.where(arr == nodata, np.nan, arr)
            src_bounds = ds.bounds
            src_crs = ds.crs
            lon_min, lat_min, lon_max, lat_max = transform_bounds(
                src_crs, "EPSG:4326",
                src_bounds.left, src_bounds.bottom, src_bounds.right, src_bounds.top,
                densify_pts=21,
            )
            # rasterio는 top-to-bottom: 첫 행이 lat_max
            bounds = (lat_max, lon_min, lat_min, lon_max)
            return arr, bounds, "local_geotiff"
    except Exception:
        return None


def _generate_ring_points(inner_xyxy: tuple[float, float, float, float],
                          outer_xyxy: tuple[float, float, float, float],
                          grid_step: float) -> np.ndarray:
    """도넛(ring) 영역을 규칙 격자로 채운 XY 좌표 반환. inner bbox 내부는 제외.

    Legacy(사각 bbox 경계): hull을 못 얻는 경우의 폴백.
    """
    ix0, iy0, ix1, iy1 = inner_xyxy
    ox0, oy0, ox1, oy1 = outer_xyxy
    xs = np.arange(ox0, ox1 + grid_step * 0.5, grid_step)
    ys = np.arange(oy0, oy1 + grid_step * 0.5, grid_step)
    gx, gy = np.meshgrid(xs, ys)
    pts = np.column_stack([gx.ravel(), gy.ravel()])
    inside = (pts[:, 0] > ix0) & (pts[:, 0] < ix1) & (pts[:, 1] > iy0) & (pts[:, 1] < iy1)
    ring_pts = pts[~inside]
    # inner 경계 자체에 정확히 놓이는 점들도 추가 (seam 정합용)
    ex = np.arange(ix0, ix1 + grid_step * 0.5, grid_step)
    edge_top = np.column_stack([ex, np.full_like(ex, iy1)])
    edge_bot = np.column_stack([ex, np.full_like(ex, iy0)])
    ey = np.arange(iy0, iy1 + grid_step * 0.5, grid_step)
    edge_left = np.column_stack([np.full_like(ey, ix0), ey])
    edge_right = np.column_stack([np.full_like(ey, ix1), ey])
    return np.vstack([ring_pts, edge_top, edge_bot, edge_left, edge_right])


def _compute_tin_hull_projected(tin_xyz_zerobased: np.ndarray,
                                origin: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """TIN 점들의 XY 컨벡스 헐을 투영 CRS 좌표로 반환.

    Returns:
        hull_xy_abs: (N,2) 반시계 순 헐 정점 XY (projected CRS, 절대좌표)
        hull_z_abs: (N,) 해당 정점의 TIN Z (절대 높이)
    """
    from scipy.spatial import ConvexHull
    xy_abs = tin_xyz_zerobased[:, :2] + origin[:2]
    hull = ConvexHull(xy_abs)
    hull_xy_abs = xy_abs[hull.vertices]
    hull_z_abs = tin_xyz_zerobased[hull.vertices, 2] + origin[2]
    return hull_xy_abs, hull_z_abs


def _densify_polygon(poly_xy: np.ndarray, max_seg_len: float) -> np.ndarray:
    """다각형 엣지를 max_seg_len 이하로 선형 보간해 세분화."""
    n = len(poly_xy)
    out = []
    for i in range(n):
        p0 = poly_xy[i]
        p1 = poly_xy[(i + 1) % n]
        seg = p1 - p0
        L = float(np.hypot(seg[0], seg[1]))
        k = max(1, int(np.ceil(L / max(max_seg_len, 1e-6))))
        out.extend(p0 + seg * t for t in np.linspace(0.0, 1.0, k, endpoint=False))
    return np.asarray(out, dtype=np.float64)


def _generate_ring_points_hull(outer_xyxy: tuple[float, float, float, float],
                               inner_hull_xy: np.ndarray,
                               grid_step: float,
                               precomputed_boundary: np.ndarray | None = None,
                               ) -> tuple[np.ndarray, int]:
    """hull 바깥쪽 격자점 + hull 엣지 세분화 점들을 결합해 링 XY 반환.

    Args:
        precomputed_boundary: (N,2) TIN mesh에 이미 존재하는 bbox 변 정점 XY.
            제공되면 densify를 건너뛰고 이 점들을 그대로 inner 경계로 사용 →
            TIN 정점과 DEM 링 정점이 **공유 XY**가 되어 T-vertex/fin 원천 제거.

    Returns:
        pts: (M, 2) 전체 링 점
        n_hull_boundary: 끝 n개가 hull 경계 점. Z 오버라이드/feather용.
    """
    ox0, oy0, ox1, oy1 = outer_xyxy
    xs = np.arange(ox0, ox1 + grid_step * 0.5, grid_step)
    ys = np.arange(oy0, oy1 + grid_step * 0.5, grid_step)
    gx, gy = np.meshgrid(xs, ys)
    pts = np.column_stack([gx.ravel(), gy.ravel()])
    # inner_hull = TIN bbox 4 모서리. bbox 내부(margin 포함)의 격자점 제거.
    # 여유 margin = grid_step*0.5 — bbox 선에 너무 가까운 격자점도 제거해 링 삼각형이
    # bbox 내부로 밀고 들어가는 현상(error1.png 톱니 fin) 차단.
    ix0 = float(np.min(inner_hull_xy[:, 0])); ix1 = float(np.max(inner_hull_xy[:, 0]))
    iy0 = float(np.min(inner_hull_xy[:, 1])); iy1 = float(np.max(inner_hull_xy[:, 1]))
    guard = grid_step * 0.5
    near_inner = (
        (pts[:, 0] > ix0 - guard) & (pts[:, 0] < ix1 + guard)
        & (pts[:, 1] > iy0 - guard) & (pts[:, 1] < iy1 + guard)
    )
    outside_grid = pts[~near_inner]

    if precomputed_boundary is not None and len(precomputed_boundary) >= 3:
        # **경계 densify 추가** — TIN bbox 공유정점 사이 간격이 외곽 grid_step 보다
        # 크면 링 Delaunay 가 bbox 를 가로질러 큰 삼각형을 만들어 TIN 영역을 덮음.
        # 따라서 공유정점을 유지하되, 인접 공유정점 사이 gap > grid_step 구간은
        # 외곽 격자와 같은 밀도로 세분화해 삽입.
        raw = np.asarray(precomputed_boundary, dtype=np.float64)
        b_tol = max(ix1 - ix0, iy1 - iy0) * 1e-4 + 1e-3
        on_bot = np.abs(raw[:, 1] - iy0) < b_tol
        on_top = np.abs(raw[:, 1] - iy1) < b_tol
        on_lft = np.abs(raw[:, 0] - ix0) < b_tol
        on_rgt = np.abs(raw[:, 0] - ix1) < b_tol
        parts: list[np.ndarray] = []
        for mask, sort_col, fixed_col, fixed_val in [
            (on_bot, 0, 1, iy0), (on_top, 0, 1, iy1),
            (on_lft, 1, 0, ix0), (on_rgt, 1, 0, ix1),
        ]:
            if mask.sum() == 0:
                continue
            side = raw[mask]
            # 변 sort + 중복 제거
            side = side[np.argsort(side[:, sort_col])]
            _, uq = np.unique(np.round(side[:, sort_col], 3), return_index=True)
            side = side[np.sort(uq)]
            parts.append(side)
            if len(side) < 2:
                continue
            main = side[:, sort_col]
            gaps = np.diff(main)
            for k, g in enumerate(gaps):
                if g > grid_step:
                    n_add = int(np.ceil(g / grid_step)) - 1
                    if n_add < 1:
                        continue
                    mids = np.linspace(main[k], main[k + 1], n_add + 2)[1:-1]
                    add = np.zeros((len(mids), 2))
                    add[:, sort_col] = mids
                    add[:, fixed_col] = fixed_val
                    parts.append(add)
        boundary_pts = np.unique(np.round(np.vstack(parts), 3), axis=0) if parts else raw
    else:
        boundary_pts = _densify_polygon(inner_hull_xy, grid_step)

    # 최종 중복 제거 (outside_grid 점 중 boundary 와 매우 가까운 것 제거 → sliver 방지)
    if len(boundary_pts) > 0 and len(outside_grid) > 0:
        from scipy.spatial import cKDTree as _cKDT
        bt = _cKDT(boundary_pts)
        d, _ = bt.query(outside_grid, k=1)
        outside_grid = outside_grid[d > grid_step * 0.25]

    n_hull_boundary = len(boundary_pts)
    all_pts = np.vstack([outside_grid, boundary_pts])
    return all_pts, n_hull_boundary


def build_extended_terrain_ring(projected_bounds: tuple[float, float, float, float],
                                origin: np.ndarray,
                                src_crs: str,
                                buffer_m: float = 1000.0,
                                tin_xyz_zerobased: np.ndarray | None = None,
                                grid_step_m: float | None = None,
                                feather_m: float = 80.0,
                                zoom: int = 13,
                                cache_dir: str | Path = "cache/dem",
                                source: str = "auto",
                                use_hull_boundary: bool = True,
                                datum_offset_override: float | None = None,
                                elev_grid_override: np.ndarray | None = None,
                                grid_bounds_override: tuple[float, float, float, float] | None = None,
                                log_fn: Callable[[str], None] = print) -> DemExtendResult:
    """DXF bbox 외곽을 실제 DEM으로 확장한 메시 생성.

    Args:
        projected_bounds: (xmin, ymin, xmax, ymax) 투영 CRS 원본 좌표.
        origin: (ox, oy, oz) — TIN 생성 시 zero-basing에 사용한 offset.
        src_crs: 투영 CRS 문자열 (예: "EPSG:5187").
        buffer_m: 외곽으로 확장할 거리(미터).
        tin_xyz_zerobased: (N,3) zero-based TIN 정점. 제공되면 수직 datum 보정
            + feathering + (use_hull_boundary=True면) 컨벡스 헐 경계 스냅에 사용.
            None이면 DEM Z 원본을 그대로 사용하고 bbox 경계 사용.
        grid_step_m: 외곽 링 격자 간격. None이면 buffer/25 기준 자동 (최소 10m).
        feather_m: inner 경계에서 이 거리 이내는 TIN Z로 블렌드.
        zoom: AWS terrarium zoom level (13 ≈ 19m/px at mid-lat).
        cache_dir: 캐시 디렉토리.
        source: "auto" | "aws_terrain" | "local_geotiff".
        use_hull_boundary: True면 링 inner 경계를 TIN 컨벡스 헐로 스냅. 헐 정점은
            TIN Z를 그대로 써 공유 정점으로 삼아 gap/T-vertex 제거. False면 사각 bbox
            경계(legacy).
        datum_offset_override: None이 아니면 이 값을 DEM Z에서 차감(vertical datum 보정).
            TIN/내부채움과 **같은 offset**을 써야 bbox seam에 Z 단차가 없음. 호출자가
            create_tin_from_dxf의 `_dem_datum_offset`을 그대로 넘기는 것을 권장.
        elev_grid_override, grid_bounds_override: None이 아니면 DEM 타일을 **재사용**
            (TIN 생성 시 이미 받은 격자 재활용 → datum 일관성 + 속도).
        log_fn: 로그 callback.

    Returns:
        DemExtendResult: 메시 + 메타데이터.
    """
    cache_dir = Path(cache_dir)
    cache_dir.mkdir(parents=True, exist_ok=True)

    xmin, ymin, xmax, ymax = projected_bounds
    inner = (xmin, ymin, xmax, ymax)
    outer = (xmin - buffer_m, ymin - buffer_m, xmax + buffer_m, ymax + buffer_m)

    if grid_step_m is None:
        grid_step_m = max(10.0, buffer_m / 25.0)

    # --- 좌표계 진단 로그 -------------------------------------------------
    # TIN/DXF는 src_crs(예: EPSG:5187), DEM 타일은 WGS84로 받아와 Z만 샘플링.
    # 결과 mesh의 XY는 src_crs 평면 그대로이므로 TIN과 동일 CRS. 혹시 변환
    # 오차가 의심되면 왕복 변환 거리로 확인 가능.
    log_fn(f"  [DEM] 좌표계: DXF={src_crs} → WGS84로 DEM 샘플링 → 결과 mesh는 {src_crs} 평면")
    try:
        _mid_x = 0.5 * (xmin + xmax); _mid_y = 0.5 * (ymin + ymax)
        _to_wgs_t = pyproj.Transformer.from_crs(src_crs, "EPSG:4326", always_xy=True)
        _to_proj_t = pyproj.Transformer.from_crs("EPSG:4326", src_crs, always_xy=True)
        _lon, _lat = _to_wgs_t.transform(_mid_x, _mid_y)
        _rx, _ry = _to_proj_t.transform(_lon, _lat)
        _rt_err = float(((_rx - _mid_x) ** 2 + (_ry - _mid_y) ** 2) ** 0.5)
        log_fn(f"  [DEM] TIN bbox: X=[{xmin:.1f}, {xmax:.1f}] Y=[{ymin:.1f}, {ymax:.1f}] "
               f"(중심 WGS84 Lon={_lon:.5f} Lat={_lat:.5f})")
        log_fn(f"  [DEM] 좌표 왕복 변환 오차: {_rt_err:.4f}m (중심점 기준, 0에 가까워야 정상)")
    except Exception as _ce:
        log_fn(f"  [DEM] 좌표계 진단 실패: {_ce}")

    # 1) 링 점 생성 — **TIN bbox를 inner polygon으로 사용** (convex hull 대신).
    #    이전 hull 모드는 TIN pts의 bbox ≠ convex hull 이라 두 메시가 hull 외곽과
    #    bbox 사이 얇은 영역을 중복 덮거나 T-vertex/fin을 만들었다. bbox를 inner로
    #    쓰면 DEM 링 inner 경계가 TIN bbox 선과 정확히 맞닿아 중복/fin 제거.
    hull_xy_abs = None
    n_hull_boundary = 0
    boundary_mode = "bbox"
    if use_hull_boundary and tin_xyz_zerobased is not None and len(tin_xyz_zerobased) >= 3:
        try:
            tin_abs = tin_xyz_zerobased + origin  # (N,3) 절대좌표
            tbx0 = float(tin_abs[:, 0].min()); tbx1 = float(tin_abs[:, 0].max())
            tby0 = float(tin_abs[:, 1].min()); tby1 = float(tin_abs[:, 1].max())
            hull_xy_abs = np.array([
                [tbx0, tby0], [tbx1, tby0], [tbx1, tby1], [tbx0, tby1],
            ], dtype=np.float64)
            # **TIN의 실제 bbox 변 정점**만 선별해 DEM 링 inner 경계로 사용 →
            # 두 메시가 동일 XY의 공유 정점을 가지므로 T-vertex/fin 원천 제거.
            bbox_tol_tin = max(tbx1 - tbx0, tby1 - tby0) * 1e-4 + 1e-3
            on_bbox_mask = (
                (np.abs(tin_abs[:, 0] - tbx0) < bbox_tol_tin)
                | (np.abs(tin_abs[:, 0] - tbx1) < bbox_tol_tin)
                | (np.abs(tin_abs[:, 1] - tby0) < bbox_tol_tin)
                | (np.abs(tin_abs[:, 1] - tby1) < bbox_tol_tin)
            )
            tin_bbox_vertices_xy = tin_abs[on_bbox_mask, :2]
            if len(tin_bbox_vertices_xy) < 4:
                # bbox 변에 정점이 너무 적으면 폴백(densify 기본 로직)
                tin_bbox_vertices_xy = None

            ring_xy, n_hull_boundary = _generate_ring_points_hull(
                outer, hull_xy_abs, grid_step_m,
                precomputed_boundary=tin_bbox_vertices_xy,
            )
            n_shared = len(tin_bbox_vertices_xy) if tin_bbox_vertices_xy is not None else 0
            boundary_mode = (f"TIN bbox 공유정점 {n_shared}개" if n_shared > 0
                             else "TIN bbox 4모서리 densify 폴백")
        except Exception as e:
            log_fn(f"  [DEM] bbox 경계 구성 실패 ({e}) — 단순 bbox 링 폴백")
            ring_xy = _generate_ring_points(inner, outer, grid_step_m)
            hull_xy_abs = None; n_hull_boundary = 0
    else:
        ring_xy = _generate_ring_points(inner, outer, grid_step_m)
    log_fn(f"  [DEM] 링 격자점: {len(ring_xy)}개 (step={grid_step_m:.1f}m, "
           f"buffer={buffer_m:.0f}m, 경계={boundary_mode})")

    # 2) WGS84로 변환
    to_wgs = pyproj.Transformer.from_crs(src_crs, "EPSG:4326", always_xy=True)
    lons, lats = to_wgs.transform(ring_xy[:, 0], ring_xy[:, 1])
    lons = np.asarray(lons); lats = np.asarray(lats)

    # 3) DEM 소스 결정 및 fetch
    resolved_source = source
    elev_grid = None
    grid_bounds = None

    # 3-0) 호출자 override — TIN densify에서 이미 받은 격자 재사용 (datum 일관성)
    if elev_grid_override is not None and grid_bounds_override is not None:
        elev_grid = elev_grid_override
        grid_bounds = grid_bounds_override
        resolved_source = "reused_from_caller"
        log_fn(f"  [DEM] 호출자 override 격자 재사용 ({elev_grid.shape[1]}x{elev_grid.shape[0]}) "
               f"— TIN 생성 때 받은 것과 **동일 datum 보장**")

    if elev_grid is None and source in ("auto", "local_geotiff"):
        local = _try_local_geotiff(cache_dir)
        if local is not None:
            elev_grid, grid_bounds, _ = local
            resolved_source = "local_geotiff"
            log_fn(f"  [DEM] 로컬 GeoTIFF 사용: cache/dem/local.tif ({elev_grid.shape[1]}x{elev_grid.shape[0]})")

    if elev_grid is None:
        ox0, oy0, ox1, oy1 = outer
        margin = max(grid_step_m * 2, 50.0)
        corners_x = np.array([ox0 - margin, ox1 + margin, ox0 - margin, ox1 + margin])
        corners_y = np.array([oy0 - margin, oy0 - margin, oy1 + margin, oy1 + margin])
        cx_lon, cx_lat = to_wgs.transform(corners_x, corners_y)
        min_lat, max_lat = float(np.min(cx_lat)), float(np.max(cx_lat))
        min_lon, max_lon = float(np.min(cx_lon)), float(np.max(cx_lon))
        elev_grid, grid_bounds = fetch_terrarium_grid(
            min_lat, min_lon, max_lat, max_lon,
            zoom=zoom, cache_dir=cache_dir, log_fn=log_fn,
        )
        resolved_source = "aws_terrain"

    # 4) 링 점 고도 샘플링
    z_dem_raw = _sample_grid_bilinear(elev_grid, grid_bounds, lats, lons)
    if np.any(np.isnan(z_dem_raw)):
        med = float(np.nanmedian(z_dem_raw))
        z_dem_raw = np.where(np.isnan(z_dem_raw), med, z_dem_raw)

    # 4a) 전역 outlier 클립 + **NaN 즉시 채움**.
    #     terrarium 디코딩 실패·타일 엣지에서 극단값/NaN 이 링 삼각형을 수십 m
    #     솟게 만들었던 현상 방지. IQR 기반 + 절대범위(−500m~9000m) 가드.
    finite = np.isfinite(z_dem_raw)
    if finite.any():
        vals = z_dem_raw[finite]
        q1, q3 = np.percentile(vals, [5, 95])
        iqr = max(q3 - q1, 1.0)
        lo = max(-500.0, q1 - iqr * 3.0)
        hi = min(9000.0, q3 + iqr * 3.0)
        med_all = float(np.median(vals))
        z_dem_raw = np.where(finite, z_dem_raw, med_all)
        z_dem_raw = np.clip(z_dem_raw, lo, hi)
    else:
        z_dem_raw = np.zeros_like(z_dem_raw)
        log_fn("  [DEM] 경고: 샘플 전체가 NaN — 0으로 채움 (DEM 타일 fetch 점검 필요)")

    # 4b) 국소 스파이크 필터 — 각 점의 Z가 반경 R 내 이웃 median 대비 3×MAD 이상
    #     벗어나면 median으로 치환. 경계 값이 튀는 원인(인접 타일 디코딩 차이)을 직접 제거.
    try:
        radius = max(grid_step_m * 3.5, 60.0)
        tree_local = cKDTree(ring_xy)
        neigh_idxs = tree_local.query_ball_point(ring_xy, r=radius)
        cleaned = z_dem_raw.copy()
        spike_count = 0
        for i, nbrs in enumerate(neigh_idxs):
            if len(nbrs) < 5:
                continue
            neigh_z = z_dem_raw[nbrs]
            med_ = float(np.median(neigh_z))
            mad_ = float(np.median(np.abs(neigh_z - med_))) + 1e-6
            if abs(z_dem_raw[i] - med_) > 4.0 * mad_:
                cleaned[i] = med_
                spike_count += 1
        if spike_count:
            log_fn(f"  [DEM] 스파이크 {spike_count}개 제거 (radius={radius:.0f}m, 4×MAD)")
        z_dem_raw = cleaned
    except Exception as _se:
        log_fn(f"  [DEM] 스파이크 필터 경고: {_se}")

    # 5) 수직 datum 보정 + feathering
    # **override 우선** — TIN 생성 시 계산한 offset을 그대로 쓰면 bbox seam에 Z 단차 0.
    vertical_offset = 0.0
    z_final = z_dem_raw.copy()
    if datum_offset_override is not None:
        vertical_offset = float(datum_offset_override)
        z_final = z_dem_raw - vertical_offset
        log_fn(f"  [DEM] 수직 datum **override** offset={vertical_offset:+.2f}m "
               f"(TIN 생성 시 계산한 값 재사용 → bbox seam Z 단차 0)")
    if tin_xyz_zerobased is not None and len(tin_xyz_zerobased) > 0:
        tin_abs = tin_xyz_zerobased + origin
        tree = cKDTree(tin_abs[:, :2])
        ring_abs = np.column_stack([ring_xy[:, 0], ring_xy[:, 1]])
        dists, idxs = tree.query(ring_abs, k=1)
        near_mask = dists < max(feather_m * 1.5, 50.0)
        if np.any(near_mask) and datum_offset_override is None:
            tin_z_near = tin_abs[idxs[near_mask], 2]
            dem_z_near = z_dem_raw[near_mask]
            diffs = dem_z_near - tin_z_near
            vertical_offset = float(np.median(diffs))
            log_fn(f"  [DEM] 수직 datum 자동 보정: offset={vertical_offset:+.2f}m "
                   f"(비교 점 {int(near_mask.sum())}개)")
            z_final = z_dem_raw - vertical_offset
        if feather_m > 0 and np.any(near_mask):
            # smoothstep 블렌드 (선형보다 부드러운 C1 전이 — 경계 Z 단차 제거).
            # t=0(경계)에서 TIN Z, t=1(feather_m 바깥)에서 DEM Z. 사이 값은
            # 3t²−2t³ 곡선으로 양 끝에서 미분 0 → 절벽 같은 Z 점프 제거.
            raw_t = np.clip(dists / feather_m, 0.0, 1.0)
            t = raw_t * raw_t * (3.0 - 2.0 * raw_t)   # smoothstep
            tin_z_all = tin_abs[idxs, 2]
            blended = (1.0 - t) * tin_z_all + t * z_final
            z_final = np.where(dists < feather_m, blended, z_final)

    # 5b) **Outer smooth blend** — outer 경계 근처를 이웃 평균으로 부드럽게 평활.
    #     이전 "하위 10% 분위수" ramp-down 은 점마다 강제 하강량이 달라 톱니 fin 을
    #     만들었음(error1.png). 이번에는:
    #       - 이웃 **평균** Z 를 타겟으로(분위수 NO)
    #       - smoothstep 가중치 (경계 1, ramp 밖 0)
    #       - 이웃 반경 = grid_step*6 로 넓게 → fin 유발하는 국소 편차 제거
    #       - 최대 하강/상승 폭 제한 (|Δ| ≤ grid_step*0.8) → 튀는 값 차단
    try:
        ox0, oy0, ox1, oy1 = outer
        dist_to_outer = np.minimum.reduce([
            ring_xy[:, 0] - ox0,
            ox1 - ring_xy[:, 0],
            ring_xy[:, 1] - oy0,
            oy1 - ring_xy[:, 1],
        ])
        dist_to_outer = np.maximum(dist_to_outer, 0.0)
        ramp_width = max(float(grid_step_m) * 4.0, float(buffer_m) * 0.20)
        in_ramp = dist_to_outer < ramp_width
        n_ramp = int(in_ramp.sum())
        if n_ramp > 0:
            tree_outer = cKDTree(ring_xy)
            r_neigh = max(float(grid_step_m) * 6.0, 120.0)
            idxs_ramp = np.where(in_ramp)[0]
            nbrs_ramp = tree_outer.query_ball_point(ring_xy[idxs_ramp], r=r_neigh)
            z_target = z_final[idxs_ramp].copy()
            for k, nb in enumerate(nbrs_ramp):
                if len(nb) < 5:
                    continue
                z_target[k] = float(np.mean(z_final[nb]))
            # smoothstep weight: outer(t=0) → w=1, ramp 밖(t=1) → w=0
            t = np.clip(dist_to_outer[idxs_ramp] / ramp_width, 0.0, 1.0)
            u = 1.0 - t
            w = u * u * (3.0 - 2.0 * u)
            z_before = z_final[idxs_ramp].copy()
            blended = (1.0 - w) * z_before + w * z_target
            # 튀는 값 가드 — 이웃 평균과의 편차를 grid_step*0.8 이내로 제한
            cap = float(grid_step_m) * 0.8
            blended = np.clip(blended, z_before - cap, z_before + cap)
            z_final[idxs_ramp] = blended
            log_fn(
                f"  [DEM] outer smooth blend {n_ramp}개 점 "
                f"(ramp폭={ramp_width:.0f}m, r={r_neigh:.0f}m, |Δ|≤{cap:.1f}m) "
                f"— 이웃 평균 블렌드, fin/톱니 제거"
            )
    except Exception as _rd:
        log_fn(f"  [DEM] outer smooth blend 경고: {_rd}")

    # 5c) bbox 경계 densified 점 Z 설정.
    #  (1) **TIN 공유 정점** (= 실제 TIN vertex 그대로) 은 TIN Z 그대로.
    #  (2) **경계 보간 densify** 신규 점(= TIN vertex 없는 bbox 변 위 점)은 양옆
    #      TIN vertex 사이 선형 보간 Z (거리 가중 k=2) 로 설정 → TIN 경계선이
    #      그대로 연장되어 seam 에서 Z 점프 0.
    if n_hull_boundary > 0 and hull_xy_abs is not None and tin_xyz_zerobased is not None:
        n_total = len(ring_xy)
        boundary_start = n_total - n_hull_boundary
        boundary_xy = ring_xy[boundary_start:]
        tin_abs_full = tin_xyz_zerobased + origin
        tin_tree_full = cKDTree(tin_abs_full[:, :2])
        d_near, i_near = tin_tree_full.query(boundary_xy, k=2)
        # k=1 이면 거의 일치하는 TIN vertex → 그대로 사용
        # k=2 가중평균으로 densify 중간 점 Z 를 선형 보간
        exact = d_near[:, 0] < 1e-3
        z_b = np.empty(len(boundary_xy), dtype=np.float64)
        z_b[exact] = tin_abs_full[i_near[exact, 0], 2]
        if (~exact).any():
            w0 = 1.0 / np.maximum(d_near[~exact, 0], 1e-6) ** 2
            w1 = 1.0 / np.maximum(d_near[~exact, 1], 1e-6) ** 2
            z0 = tin_abs_full[i_near[~exact, 0], 2]
            z1 = tin_abs_full[i_near[~exact, 1], 2]
            z_b[~exact] = (w0 * z0 + w1 * z1) / (w0 + w1)
        z_final[boundary_start:] = z_b

    # 5d) **최종 NaN/무한 가드** — 어떤 경로로 들어온 NaN도 여기서 확정적으로 제거.
    #     인접 이웃의 median 으로 채움. 이웃 없으면 전체 median.
    bad = ~np.isfinite(z_final)
    if bad.any():
        good_idx = np.where(~bad)[0]
        if len(good_idx) > 0:
            tree_nan = cKDTree(ring_xy[good_idx])
            _, j = tree_nan.query(ring_xy[bad], k=1)
            z_final[bad] = z_final[good_idx[j]]
        else:
            z_final[bad] = 0.0
        log_fn(f"  [DEM] NaN/Inf {int(bad.sum())}개 이웃 값으로 채움")

    # 5e) **링 Z Laplacian 1pass** — 톱니(fin) 마지막 평활. TIN 공유 경계 정점은
    #     고정(Laplacian 에서 제외) → TIN–링 seam Z 는 변경되지 않음.
    try:
        tree_ring = cKDTree(ring_xy)
        r_sm = max(float(grid_step_m) * 2.2, 40.0)
        fixed_mask = np.zeros(len(ring_xy), dtype=bool)
        if n_hull_boundary > 0:
            fixed_mask[len(ring_xy) - n_hull_boundary:] = True
        movable = np.where(~fixed_mask)[0]
        nbrs = tree_ring.query_ball_point(ring_xy[movable], r=r_sm)
        z_new = z_final.copy()
        for k, nb in enumerate(nbrs):
            if len(nb) < 4:
                continue
            vid = movable[k]
            nb_arr = np.asarray(nb, dtype=np.int64)
            nb_arr = nb_arr[nb_arr != vid]
            if len(nb_arr) == 0:
                continue
            z_new[vid] = 0.5 * z_final[vid] + 0.5 * float(np.mean(z_final[nb_arr]))
        z_final = z_new
    except Exception as _ls:
        log_fn(f"  [DEM] 링 Laplacian 평활 경고: {_ls}")

    # 6) Zero-basing
    ring_xyz = np.column_stack([
        ring_xy[:, 0] - origin[0],
        ring_xy[:, 1] - origin[1],
        z_final - origin[2],
    ])

    # 7) Delaunay (XY 기준) → 도넛 triangulation
    try:
        tri = Delaunay(ring_xyz[:, :2])
    except Exception as e:
        raise RuntimeError(f"링 메시 Delaunay 실패: {e}") from e

    # inner(TIN bbox) 내부 삼각형 제거 — **3중 가드**로 침범 원천 차단.
    #   (a) centroid 가 inner bbox 내부
    #   (b) **세 정점 중 하나라도 inner bbox strict 내부** (≥ strict_tol 안쪽)
    #   (c) **엣지 중점** 셋 중 하나라도 inner bbox strict 내부
    # 이전에는 (a)+(b) 만으로 bbox 선을 따라 "세 정점 전부 경계 위"인 납작한
    # 삼각형이 남아 엣지 중점이 bbox 내부로 진입하며 TIN 위를 덮었음(error.png
    # 침범 증상). (c) 추가로 완전 차단.
    cx0 = inner[0] - origin[0]; cy0 = inner[1] - origin[1]
    cx1 = inner[2] - origin[0]; cy1 = inner[3] - origin[1]
    bbox_w = max(cx1 - cx0, cy1 - cy0)
    strict_tol = bbox_w * 1e-4 + 1e-3  # 수치 오차 가드 (≈mm 단위)

    v0 = ring_xyz[tri.simplices[:, 0], :2]
    v1 = ring_xyz[tri.simplices[:, 1], :2]
    v2 = ring_xyz[tri.simplices[:, 2], :2]
    centroids = (v0 + v1 + v2) / 3.0
    mid01 = (v0 + v1) * 0.5
    mid12 = (v1 + v2) * 0.5
    mid20 = (v2 + v0) * 0.5

    def _strict_in_bbox(p: np.ndarray) -> np.ndarray:
        return (
            (p[:, 0] > cx0 + strict_tol) & (p[:, 0] < cx1 - strict_tol)
            & (p[:, 1] > cy0 + strict_tol) & (p[:, 1] < cy1 - strict_tol)
        )

    # (a) centroid
    inside_inner_centroid = _strict_in_bbox(centroids - np.array([[0, 0]]))  # strict only
    # centroid 는 strict 내부까지 너무 보수적이면 seam 삼각형도 지울 수 있으므로
    # 너그럽게: bbox 완전 내부만.
    inside_inner_centroid = (
        (centroids[:, 0] > cx0 + strict_tol) & (centroids[:, 0] < cx1 - strict_tol)
        & (centroids[:, 1] > cy0 + strict_tol) & (centroids[:, 1] < cy1 - strict_tol)
    )

    # (b) 세 정점 중 하나라도 strict 내부
    vx = ring_xyz[:, 0]; vy = ring_xyz[:, 1]
    strict_inner_vert = (
        (vx > cx0 + strict_tol) & (vx < cx1 - strict_tol)
        & (vy > cy0 + strict_tol) & (vy < cy1 - strict_tol)
    )
    has_strict_vertex = (
        strict_inner_vert[tri.simplices[:, 0]]
        | strict_inner_vert[tri.simplices[:, 1]]
        | strict_inner_vert[tri.simplices[:, 2]]
    )

    # (c) 엣지 중점 중 하나라도 strict 내부
    edge_mid_inside = (
        _strict_in_bbox(mid01) | _strict_in_bbox(mid12) | _strict_in_bbox(mid20)
    )

    inside_inner = inside_inner_centroid | has_strict_vertex | edge_mid_inside
    n_cent = int(inside_inner_centroid.sum())
    n_strict = int(has_strict_vertex.sum())
    n_emid = int(edge_mid_inside.sum())
    tris_keep = tri.simplices[~inside_inner]
    log_fn(f"  [DEM] inner 제거 삼각형: centroid {n_cent}개 + strict-vertex {n_strict}개 "
           f"+ edge-midpoint {n_emid}개 (tol={strict_tol:.3f}m) — TIN 침범 3중 차단")

    # 링 삼각형 벽 컷 — **2단계 보호 + 추가 완화**.
    # 이전: slope_ratio>2.5 & Z스팬>10m 도 실제 72° 산사면/급한 골짜기를 자르면서
    #       접합점에 삼각 구멍(= error.png의 두 번째 이슈)을 냈다.
    # 이번:
    #   (1) **TIN 공유 경계 정점을 가진 삼각형은 컷 제외** — bbox seam 삼각형은
    #       DEM 급경사 한 번에 날아가면 보이는 **접합점 구멍**이 된다. 보호.
    #   (2) 기준 상향: slope_ratio > 4.0 (≈76°) AND Z스팬 > 30m — 자연 급사면보다
    #       가파르고, 30m 이상 Z 점프가 있어야만 진짜 "수직 벽" 후보로 간주.
    #   (3) 추가 가드: 에지 길이 e_max가 너무 짧지 않아야 함(<5m는 제외: 극미세
    #       sliver는 별도 수치 오차, 벽 아님).
    if len(tris_keep) > 0:
        zs = ring_xyz[tris_keep][:, :, 2]
        z_span_tri = zs.max(axis=1) - zs.min(axis=1)
        pxy = ring_xyz[tris_keep][:, :, :2]
        e01 = np.linalg.norm(pxy[:, 0] - pxy[:, 1], axis=1)
        e12 = np.linalg.norm(pxy[:, 1] - pxy[:, 2], axis=1)
        e20 = np.linalg.norm(pxy[:, 2] - pxy[:, 0], axis=1)
        e_max_tri = np.maximum(np.maximum(e01, e12), e20)
        e_max_safe = np.maximum(e_max_tri, 1e-6)
        slope_ratio_tri = z_span_tri / e_max_safe

        # (1) seam 보호 마스크 — TIN bbox 변 공유 정점(= boundary_pts, 마지막
        #     n_hull_boundary 개)을 가진 삼각형은 "접합부"라 컷 제외.
        protect_vertex = np.zeros(len(ring_xyz), dtype=bool)
        if n_hull_boundary > 0:
            protect_vertex[len(ring_xyz) - n_hull_boundary:] = True
        pv0 = protect_vertex[tris_keep[:, 0]]
        pv1 = protect_vertex[tris_keep[:, 1]]
        pv2 = protect_vertex[tris_keep[:, 2]]
        seam_tri = pv0 | pv1 | pv2

        # (2+3) 완화된 기준
        wall_mask = (slope_ratio_tri > 4.0) & (z_span_tri > 30.0) & (e_max_tri > 5.0) & (~seam_tri)
        if wall_mask.any():
            tris_keep = tris_keep[~wall_mask]
            log_fn(f"  [DEM] 링 벽 삼각형 {int(wall_mask.sum())}개 제거 "
                   f"(slope_ratio>4.0(≈76°) & Z스팬>30m, seam {int(seam_tri.sum())}개 보호) "
                   f"— 자연 급사면 보존, 접합점 구멍 방지")
        else:
            log_fn(f"  [DEM] 링 벽 컷: 대상 없음 (seam {int(seam_tri.sum())}개 보호) "
                   f"— 자연 경사면 유지")

    if len(tris_keep) == 0:
        raise RuntimeError("링 삼각형이 하나도 남지 않았습니다 (buffer_m/grid_step_m 확인).")

    faces = np.column_stack([np.full(len(tris_keep), 3), tris_keep])
    mesh = pv.PolyData(ring_xyz, faces)
    mesh["Elevation"] = ring_xyz[:, 2]

    info = (f"source={resolved_source}, buffer={buffer_m:.0f}m, step={grid_step_m:.1f}m, "
            f"feather={feather_m:.0f}m, offset={vertical_offset:+.2f}m, "
            f"boundary={boundary_mode}, pts={len(ring_xyz)}, tris={len(tris_keep)}")
    log_fn(f"  [DEM] 링 메시 완료: {info}")

    return DemExtendResult(
        mesh=mesh,
        n_points=len(ring_xyz),
        n_faces=len(tris_keep),
        buffer_m=buffer_m,
        source=resolved_source,
        vertical_offset_m=vertical_offset,
        zoom=zoom,
        grid_step_m=grid_step_m,
        feather_m=feather_m,
        info=info,
    )


if __name__ == "__main__":
    # 스모크 테스트: 사연댐 대략 좌표(WGS84) 주변
    # 투영 좌표로 변환해서 테스트
    to_proj = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:5187", always_xy=True)
    cx, cy = to_proj.transform(129.10, 35.54)  # 대략 사연댐 부근
    half = 1000.0
    bounds = (cx - half, cy - half, cx + half, cy + half)
    origin = np.array([cx - half, cy - half, 100.0])
    result = build_extended_terrain_ring(
        projected_bounds=bounds,
        origin=origin,
        src_crs="EPSG:5187",
        buffer_m=500.0,
        grid_step_m=50.0,
        feather_m=0.0,
        log_fn=print,
    )
    print(result.info)
