Aislo/resources/tester/test_b04_contour_equivalence.py

"""실데이터 등가성 — 래스터화가 기존 griddata 경로와 같은 표고를 내는지.

PLAN 2026-08-17 완료 조건 검증.
- `meshfree`: Delaunay를 그대로 쓰므로 표고 완전 일치
- `tin`: 저장 faces와 재Delaunay는 다른 삼각망이라 차이가 남는다. 그 차이가
  **보간 버그가 아니라 삼각망 선택 차이**임을 대조로 분리한다.

**주의: griddata가 호출당 3분 걸려 이 파일 전체 실행에 약 8분이 든다.**
실제 프로젝트 산출물이 없으면 skip.
"""

import time
from pathlib import Path

import numpy as np
import pytest
from scipy.interpolate import griddata
from scipy.spatial import Delaunay

from B04_PreProcess.B04_PreProcess_Engine_Contour import (
    _grid_axes,
    extract_contours,
    rasterize_triangle_mesh,
)

MODELS = Path(
    "C:/Program_coding/임도설계 및 견적자동화 프로그램 개발/storage/1/3/"
    "f9543035-4a91-4df0-b467-cdc515723507/B04_PreProcess/models"
)


def _load(name):
    path = MODELS / name
    if not path.exists():
        pytest.skip(f"실데이터 없음: {path}")
    return np.load(path)


def _axes_for(xy):
    return _grid_axes(
        float(np.min(xy[:, 0])),
        float(np.max(xy[:, 0])),
        float(np.min(xy[:, 1])),
        float(np.max(xy[:, 1])),
        1.0,
    )


def _compare(old, new):
    return {
        "only_old": int((np.isfinite(old) & ~np.isfinite(new)).sum()),
        "only_new": int((np.isfinite(new) & ~np.isfinite(old)).sum()),
        "both": np.isfinite(old) & np.isfinite(new),
    }


def test_meshfree_matches_griddata_exactly():
    """meshfree — 같은 Delaunay를 쓰므로 표고가 완전히 일치해야 한다."""
    points = _load("meshfree_csf.npz")["points"]
    x_coords, y_coords = _axes_for(points)
    xx, yy = np.meshgrid(x_coords, y_coords)

    triangulation = Delaunay(np.asarray(points[:, :2], dtype=np.float64))
    new = rasterize_triangle_mesh(points, triangulation.simplices, x_coords, y_coords)
    old = griddata(points[:, :2], points[:, 2], (xx, yy), method="linear")

    result = _compare(old, new)
    assert result["both"].sum() > 10_000, "비교 대상 셀이 너무 적다"
    assert result["only_old"] == 0, "기존이 채우던 셀이 비었다"
    assert result["only_new"] == 0, "없던 셀이 생겼다"
    assert np.abs(old[result["both"]] - new[result["both"]]).max() < 1e-9


def test_rasterize_matches_griddata_on_same_triangulation():
    """래스터화 정확성 — 같은 재Delaunay를 넣으면 griddata와 완전 일치.

    이 대조가 통과하면 tin에서 남는 차이는 전부 삼각망 선택 차이다.
    """
    vertices = _load("tin_csf.npz")["vertices"]
    x_coords, y_coords = _axes_for(vertices)
    xx, yy = np.meshgrid(x_coords, y_coords)

    triangulation = Delaunay(np.asarray(vertices[:, :2], dtype=np.float64))
    control = rasterize_triangle_mesh(vertices, triangulation.simplices, x_coords, y_coords)
    old = griddata(vertices[:, :2], vertices[:, 2], (xx, yy), method="linear")

    result = _compare(old, control)
    assert result["only_old"] == 0 and result["only_new"] == 0
    assert np.abs(old[result["both"]] - control[result["both"]]).max() < 1e-9


def test_tin_difference_comes_only_from_triangulation():
    """tin 저장 faces — griddata 대비 차이가 재Delaunay 대비 차이와 같아야 한다.

    두 대조가 일치하면 차이의 원인은 삼각망 하나뿐이다. 저장 faces는 긴 변 제거·
    외곽 클리핑이 반영된 실제 TIN이므로(`ModelBuild.py:38-51`) 재Delaunay보다 덮는
    면적이 작고, 그 몫이 `only_old`로 나온다.
    """
    data = _load("tin_csf.npz")
    vertices, faces = data["vertices"], data["faces"]
    x_coords, y_coords = _axes_for(vertices)
    xx, yy = np.meshgrid(x_coords, y_coords)

    stored = rasterize_triangle_mesh(vertices, faces, x_coords, y_coords)
    triangulation = Delaunay(np.asarray(vertices[:, :2], dtype=np.float64))
    control = rasterize_triangle_mesh(vertices, triangulation.simplices, x_coords, y_coords)
    old = griddata(vertices[:, :2], vertices[:, 2], (xx, yy), method="linear")

    vs_griddata = _compare(old, stored)
    vs_control = _compare(control, stored)
    assert vs_griddata["only_old"] == vs_control["only_old"]
    assert vs_griddata["only_new"] == vs_control["only_new"] == 0

    diff = np.abs(old[vs_griddata["both"]] - stored[vs_griddata["both"]])
    changed = diff > 1e-6
    print(
        f"\n[tin] 공통셀={int(vs_griddata['both'].sum()):,}  "
        f"저장 faces가 안 덮는 셀={vs_griddata['only_old']:,}\n"
        f"      값 다른 셀={int(changed.sum()):,} ({changed.mean():.2%})  "
        f"최대={diff.max():.4f}m  차이셀 평균={diff[changed].mean():.4f}m"
    )
    # 대각선 선택이 다른 사각형에서만 어긋나므로 국소 지형 기복을 넘지 않아야 한다.
    assert diff.max() < 1.0, "표고차 1m 초과 — 삼각망 차이로 설명 불가"


@pytest.mark.parametrize(
    "npz_name,representation,cached_name,exact",
    [
        ("tin_csf.npz", "triangular_mesh", "contour_csf_tin_1.0m.json", False),
        ("meshfree_csf.npz", "meshfree_surfels", "contour_csf_meshfree_1.0m.json", True),
        ("dtm_csf.npz", "regular_grid", "contour_csf_dtm_1.0m.json", True),
    ],
)
def test_contour_output_matches_cached(npz_name, representation, cached_name, exact):
    """최종 산출물 대조 — 캐시된 등고선 JSON과 레벨·정점 수를 맞춰 본다."""
    import json

    npz_path = MODELS / npz_name
    cached_path = MODELS / cached_name
    if not npz_path.exists() or not cached_path.exists():
        pytest.skip("실데이터 없음")

    started = time.perf_counter()
    lines = extract_contours(npz_path, representation, 1.0)
    elapsed = time.perf_counter() - started
    with open(cached_path, encoding="utf-8") as handle:
        cached = json.load(handle)["contours"]

    new_levels = {round(float(item["level"]), 3) for item in lines}
    old_levels = {round(float(item["level"]), 3) for item in cached}
    new_points = sum(len(item["coordinates"]) for item in lines)
    old_points = sum(len(item["coordinates"]) for item in cached)
    print(
        f"\n[{npz_name}] {elapsed:.2f}초  세그먼트 {len(lines):,}/{len(cached):,}  "
        f"정점 {new_points:,}/{old_points:,}"
    )

    assert new_levels == old_levels, "등고선 레벨 집합이 바뀌었다"
    if exact:
        assert len(lines) == len(cached)
        assert new_points == old_points
    else:
        # tin은 삼각망 차이만큼 경계 세그먼트가 미세하게 달라진다.
        assert abs(new_points - old_points) / old_points < 0.01


def test_contour_extraction_is_fast():
    """가속이 목적이므로 속도를 회귀로 못박는다 (변경 전 tin 119초·meshfree 302초)."""
    for npz_name, representation in (
        ("tin_csf.npz", "triangular_mesh"),
        ("meshfree_csf.npz", "meshfree_surfels"),
    ):
        path = MODELS / npz_name
        if not path.exists():
            pytest.skip("실데이터 없음")
        started = time.perf_counter()
        extract_contours(path, representation, 1.0)
        elapsed = time.perf_counter() - started
        print(f"\n[속도] {npz_name} {elapsed:.2f}초")
        assert elapsed < 30.0, f"{npz_name} {elapsed:.1f}초 — 가속 실패"