s-canvas/structure_placement.py

"""구조물 메쉬를 지형 위에 배치하는 유틸.

구조물 빌더들(intake_tower, valve_chamber, retaining_wall 등)은 구조물을
**원점 중심 + Z=구조물 설계 EL** 로컬 좌표계에서 생성한다.

이를 전체계획 평면도의 해당 구조물 위치에 배치하려면:
  1. XY: 평면도 centroid로 평행이동
  2. 회전: 평면 오리엔테이션(예: 옹벽 길이축)으로 Z축 기준 회전
  3. Z: 옵션1) 구조물 bottom_el 그대로 유지 (설계 EL 기준)
        옵션2) TIN 표면 고도로 이동 (지형에 내려앉힘)

사용법:
    from structure_placement import apply_placement

    placed = apply_placement(
        meshes=template.build_meshes(params),
        plan_centroid=(x, y),
        rotation_deg=0.0,
        z_mode="design",  # or "terrain"
        terrain_mesh=tin_mesh,  # z_mode="terrain" 시 필요
        z_offset=0.0,
    )
"""

from __future__ import annotations

import math

import numpy as np
import pyvista as pv


def rotate_points_z(points: np.ndarray, angle_rad: float,
                     cx: float = 0, cy: float = 0) -> np.ndarray:
    """Z축 기준 회전 (각 points[i] = [x, y, z])."""
    cos_a = math.cos(angle_rad)
    sin_a = math.sin(angle_rad)
    out = points.copy()
    dx = out[:, 0] - cx
    dy = out[:, 1] - cy
    out[:, 0] = cx + dx * cos_a - dy * sin_a
    out[:, 1] = cy + dx * sin_a + dy * cos_a
    return out


def interpolate_terrain_z(tin_mesh: pv.PolyData, xy: tuple[float, float],
                           origin: np.ndarray = None) -> float | None:
    """TIN 표면의 (x, y) 위치에서 Z값 보간.

    Args:
        tin_mesh: PyVista PolyData (TIN 메쉬, 로컬 좌표)
        xy: 월드 좌표 (x, y)
        origin: TIN의 원점 보정 offset [x, y, z] (scanvas_maker.self.origin)

    Returns:
        보간된 Z값 (월드 EL, m) 또는 None (범위 밖)
    """
    if tin_mesh is None:
        return None

    try:
        from scipy.interpolate import LinearNDInterpolator
        pts = np.array(tin_mesh.points)

        # origin 보정: xy를 TIN 로컬 좌표계로 변환
        if origin is not None:
            local_x = xy[0] - origin[0]
            local_y = xy[1] - origin[1]
        else:
            local_x, local_y = xy[0], xy[1]

        interp = LinearNDInterpolator(pts[:, :2], pts[:, 2])
        z_local = interp(local_x, local_y)

        if np.isnan(z_local):
            # 범위 밖 → 중앙값 사용
            z_local = float(np.median(pts[:, 2]))

        # origin Z 복원 (월드 EL로)
        if origin is not None:
            return float(z_local) + origin[2]
        return float(z_local)
    except Exception:
        return None


def apply_placement(
    meshes: list[tuple[pv.PolyData, str, float]],
    plan_centroid: tuple[float, float],
    rotation_deg: float = 0.0,
    z_mode: str = "design",
    terrain_mesh: pv.PolyData | None = None,
    terrain_origin: np.ndarray | None = None,
    z_offset: float = 0.0,
    structure_bottom_el: float | None = None,
    skip_ground: bool = False,
    scale: float = 1.0,
    skip_terrain: bool = False,
    pad_surface_z: float | None = None,
    embed_offset: float = 0.02,
) -> list[tuple[pv.PolyData, str, float]]:
    """구조물 메쉬 리스트를 지형 위 해당 위치에 배치.

    Args:
        meshes: 템플릿이 생성한 로컬 좌표 메쉬들
        plan_centroid: 평면도에서 구조물의 중심 월드 좌표 (x, y)
        rotation_deg: Z축 기준 회전 (도). 옹벽 등 선형 구조물 방향 맞춤용
        z_mode:
          - "design": 구조물 설계 EL 그대로 (Z 이동 없음)
          - "terrain": TIN 표면 고도에 bottom_el이 오도록 이동
          - "offset": z_offset만큼 이동
        terrain_mesh: z_mode="terrain"일 때 사용할 TIN
        terrain_origin: TIN의 origin 보정값
        z_offset: z_mode="offset"일 때 적용할 Z 이동량
        structure_bottom_el: 구조물 바닥 EL (z_mode="terrain"에서 지형에 맞출 때)
        skip_ground: True면 구조물 자체 "ground" 메쉬 제외 (지형 TIN과 중복 방지)
        scale: XY 방향 균등 스케일 (Geo-Referencing 결과). 1.0이면 무변화.
               Z는 의도적으로 유지 (구조물 설계 EL 보존).
        skip_terrain: True면 물/지면/에이프런/뒤채움 관련 메쉬 전부 제외
                      (geo_referencing.EXCLUDE_COLORS 집합 사용)

    Returns:
        변환된 (mesh, color, opacity) 리스트
    """
    if not meshes:
        return []

    # 모든 좌표를 TIN 로컬 좌표계로 통일
    # 템플릿 메쉬: XY=(0,0) 중심, Z=설계 EL (월드)
    # TIN 메쉬: origin이 빠진 로컬 좌표
    # → plan_centroid(월드) - origin = 로컬 배치 위치
    origin = terrain_origin if terrain_origin is not None else np.zeros(3)

    # XY 이동량: 월드 → 로컬
    dx = plan_centroid[0] - origin[0]
    dy = plan_centroid[1] - origin[1]

    # Z 이동량
    dz = 0.0
    if z_mode == "terrain" and structure_bottom_el is not None:
        if pad_surface_z is not None:
            z_local = float(pad_surface_z)
        elif terrain_mesh is not None:
            # TIN 로컬 좌표에서 직접 보간 (interpolate_terrain_z는 월드 EL 반환 → 사용 안 함)
            from scipy.interpolate import LinearNDInterpolator
            tin_pts = np.array(terrain_mesh.points)
            interp = LinearNDInterpolator(tin_pts[:, :2], tin_pts[:, 2])
            z_local = float(interp(dx, dy))
            if np.isnan(z_local):
                z_local = float(np.median(tin_pts[:, 2]))
        else:
            z_local = None
        if z_local is not None:
            # 구조물 바닥 EL(설계 월드값)을 TIN 로컬 표면에 맞춤
            # + embed_offset만큼 아래로 → TIN이 구조물 underside를 가림
            dz = z_local - structure_bottom_el - float(embed_offset)
    elif z_mode == "offset":
        dz = z_offset

    angle_rad = math.radians(rotation_deg)

    out = []
    ground_colors = {"#8B7D6B", "#7F6F5F"}  # ground 색상

    # 물/지면류 전체 제외용 색상 집합 (geo_referencing.EXCLUDE_COLORS와 동일)
    terrain_exclude_colors = {
        "#3A7AA8", "#7F6F5F", "#8B7D6B", "#9A968C", "#8B7355",
    }

    def _norm_color(c):
        if not isinstance(c, str):
            return ""
        s = c.strip().upper()
        if not s.startswith("#"):
            s = "#" + s
        return s

    # CW quad 보정 — apply_placement은 centroid/rotation 기반 폴백이므로
    # plan_centroid만 주어진 경우엔 Y-flip이 무의미함. 하지만 rotation_deg가
    # CW convention으로 주어졌다고 가정하면 y_sign 적용.
    # 단순화를 위해 apply_placement은 현재 CW Y-flip은 적용하지 않음.

    for mesh, color, opacity in meshes:
        norm = _norm_color(color)
        # ground 메쉬 건너뛰기 (지형이 있으면 중복)
        if skip_ground and norm in {c.upper() for c in ground_colors}:
            continue
        # 물/지면/apron/backfill 전부 건너뛰기 (미리보기·최종 배치용)
        if skip_terrain and norm in {c.upper() for c in terrain_exclude_colors}:
            continue

        try:
            new_pts = np.array(mesh.points).copy()

            # 1) XY 스케일 (Z 보존)
            if abs(scale - 1.0) > 1e-6:
                new_pts[:, 0] *= scale
                new_pts[:, 1] *= scale

            # 2) 회전 (원점 기준)
            if abs(rotation_deg) > 0.01:
                new_pts = rotate_points_z(new_pts, angle_rad, cx=0, cy=0)

            # 3) 평행이동 (XY + Z)
            new_pts[:, 0] += dx
            new_pts[:, 1] += dy
            new_pts[:, 2] += dz

            new_mesh = mesh.copy()
            new_mesh.points = new_pts
            out.append((new_mesh, color, opacity))
        except Exception:
            # 개별 메쉬 변환 실패시 건너뜀
            continue

    return out


def fit_meshes_to_quad(
    meshes: list[tuple[pv.PolyData, str, float]],
    quad_world_pts: list,
    terrain_mesh: pv.PolyData | None = None,
    terrain_origin: np.ndarray | None = None,
    structure_bottom_el: float | None = None,
    z_mode: str = "terrain",
    z_offset: float = 0.0,
    skip_ground: bool = True,
    skip_terrain: bool = True,
    scale_mode: str = "none",
    pad_surface_z: float | None = None,
    embed_offset: float = 0.02,
    detail_quad_pts: list | None = None,
    plan_frame_angle_deg: float = 0.0,
    flip_y_for_cw_quad: bool = True,
) -> list[tuple[pv.PolyData, str, float]]:
    """3D 메쉬를 TIN 평면도의 4점 사각형에 **강제 맞춤** (anisotropic scale).

    동작(모든 모드 공통):
      1. 메쉬 XY bbox 중심을 원점으로 이동
      2. scale_mode에 따라 스케일 적용
      3. quad 첫 엣지 각도로 XY 회전
      4. quad 중심(TIN 로컬)으로 XY 이동
      5. z_mode에 따라 Z 이동 (기본: 메쉬 z_min을 TIN 표면에 맞춤)

    Args:
        meshes: 템플릿 빌드 결과 (mesh, color_hex, opacity)
        quad_world_pts: TIN 평면도 4점 (시계방향, 월드 좌표)
        terrain_mesh: TIN PolyData (z_mode="terrain"일 때 표면 Z 보간)
        terrain_origin: TIN 로컬 좌표 변환용 origin (self.origin)
        structure_bottom_el: 호환용 (실제로는 스케일 후 mesh z_min 사용)
        z_mode: "terrain"(표면 안착) | "design"(Z 유지) | "offset"(z_offset)
        skip_ground: True면 ground 색상 메쉬 제외
        skip_terrain: True면 물/지면/apron/backfill 색상 전체 제외
        scale_mode:
          - "none" (기본/권장): 스케일 안 함. 구조물 설계 크기 유지.
          - "xy_only": quad에 맞춰 X/Y만 anisotropic 스케일 (Z 보존).
          - "xyz_uniform": sqrt(sx·sy) 균등 스케일을 X/Y/Z 모두에 적용
            + [0.01, 100] 범위로 clamp (extreme 단위 mismatch 방지)
        pad_surface_z: 굴착 pad 표면 Z (TIN 로컬). 지정 시 quad 중심 보간값 대신
            이 값을 z_surface로 사용 (굴착 영역 내부는 평탄 pad이므로 quad 중심
            보간이 아주 약간 흔들릴 수 있는 경우에도 안정).
        embed_offset: 구조물 바닥 z_min을 표면보다 이만큼 **아래로** 내림 (m).
            TIN(pad 표면)이 구조물의 아래면보다 살짝 위에 있어야 아래에서 위로
            봤을 때 TIN이 underside를 가림 (z-fighting 및 관통 방지).
            기본 0.02m(2cm) — bedding/거푸집 두께 수준으로 물리적으로도 자연스러움.
        detail_quad_pts: 구조물 상세 평면도의 4점 (사용자가 시계방향으로 picked).
            지정 시 detail ↔ TIN의 **상대 회전**을 사용하여 detail이 DXF 원본에서
            비수평으로 그려진 경우도 올바르게 처리. 없으면 절대 angle(TIN only) 폴백.
        plan_frame_angle_deg: detail DXF 내에서 mesh의 +X축이 향하는 각도(도).
            파서가 plan outline의 PCA 주축을 추출해 전달. mesh는 이 각도만큼 먼저
            회전되어 detail의 span 방향에 정렬된 후 detail→TIN 전체 변환이 적용됨.
            0 = 이미 detail의 +X를 따르는 경우.
        flip_y_for_cw_quad: True면 회전 전 mesh의 Y를 반전. 사용자가 picks를
            **시계방향**으로 찍었을 때 기본 회전만으로는 mesh의 +Y가 quad의 edge 3→0
            방향(반시계)으로 가버려 앞뒤가 뒤집힌다. Y-flip으로 mesh +Y가 edge 1→2
            방향(사용자 의도의 "downward/downstream")을 따르도록 보정.

    Returns:
        변환된 (mesh, color, opacity) 리스트 (좌표 = TIN 로컬)
    """
    if not meshes or not quad_world_pts or len(quad_world_pts) < 4:
        return []

    # 색상 필터
    ground_colors = {"#8B7D6B", "#7F6F5F"}
    terrain_colors = {"#3A7AA8", "#7F6F5F", "#8B7D6B", "#9A968C", "#8B7355"}
    exclude = set()
    if skip_ground:
        exclude |= ground_colors
    if skip_terrain:
        exclude |= terrain_colors
    exclude_upper = {c.upper() for c in exclude}

    def _norm(c):
        if not isinstance(c, str):
            return ""
        s = c.strip().upper()
        if not s.startswith("#"):
            s = "#" + s
        return s

    filtered = [(m, c, o) for (m, c, o) in meshes
                if _norm(c) not in exclude_upper]
    if not filtered:
        return []

    # 전체 메쉬 XY/Z bounding box (모든 컴포넌트 통합)
    all_pts_list = []
    for m, _, _ in filtered:
        try:
            pts = np.asarray(m.points, dtype=np.float64)
            if pts.size:
                all_pts_list.append(pts)
        except Exception:
            continue
    if not all_pts_list:
        return []
    all_pts = np.concatenate(all_pts_list, axis=0)
    mesh_xmin, mesh_ymin = float(all_pts[:, 0].min()), float(all_pts[:, 1].min())
    mesh_xmax, mesh_ymax = float(all_pts[:, 0].max()), float(all_pts[:, 1].max())
    aggregate_z_min = float(all_pts[:, 2].min())  # 모든 컴포넌트 중 최저 Z (상대 Z 보존용)
    mesh_cx = (mesh_xmin + mesh_xmax) / 2.0
    mesh_cy = (mesh_ymin + mesh_ymax) / 2.0
    mesh_w = mesh_xmax - mesh_xmin
    mesh_h = mesh_ymax - mesh_ymin

    if mesh_w < 1e-6 or mesh_h < 1e-6:
        return filtered

    # TIN quad (월드)
    quad = np.asarray(quad_world_pts[:4], dtype=np.float64)
    q_center_world = quad.mean(axis=0)

    # 월드 → TIN 로컬 (terrain_origin 차감)
    origin = np.asarray(terrain_origin if terrain_origin is not None else np.zeros(3),
                         dtype=np.float64)
    q_center_local_x = float(q_center_world[0] - origin[0])
    q_center_local_y = float(q_center_world[1] - origin[1])

    edge_01 = quad[1] - quad[0]
    edge_12 = quad[2] - quad[1]
    q_w = float(np.linalg.norm(edge_01))
    q_h = float(np.linalg.norm(edge_12))
    if q_w < 1e-6 or q_h < 1e-6:
        return filtered

    # TIN quad 첫 엣지 각도 (반시계 양수)
    tin_angle = math.atan2(float(edge_01[1]), float(edge_01[0]))

    # Detail quad의 첫 엣지 각도 (있으면 상대회전 계산용)
    if detail_quad_pts and len(detail_quad_pts) >= 4:
        d_quad = np.asarray(detail_quad_pts[:4], dtype=np.float64)
        d_edge_01 = d_quad[1] - d_quad[0]
        if np.linalg.norm(d_edge_01) >= 1e-6:
            detail_angle = math.atan2(float(d_edge_01[1]), float(d_edge_01[0]))
        else:
            detail_angle = 0.0
    else:
        # 폴백: detail의 +X가 edge 0→1과 같다고 가정 (기존 동작)
        detail_angle = 0.0

    # mesh가 detail 프레임 내에서 frame_angle만큼 회전되어 있으면 그만큼 상쇄
    # 최종 회전 = (TIN 각도) - (Detail 각도) - (mesh의 detail 내 회전)
    q_angle = tin_angle - detail_angle - math.radians(plan_frame_angle_deg)
    cos_a = math.cos(q_angle)
    sin_a = math.sin(q_angle)
    # CW quad 보정: True면 mesh Y 반전 후 회전 (mesh +Y가 edge 1→2 따라가도록)
    y_sign = -1.0 if flip_y_for_cw_quad else 1.0

    # Scale 선택
    raw_sx = q_w / mesh_w
    raw_sy = q_h / mesh_h
    if scale_mode == "xy_only":
        scale_x, scale_y, scale_z_val = raw_sx, raw_sy, 1.0
    elif scale_mode == "xyz_uniform":
        s = math.sqrt(abs(raw_sx * raw_sy))
        s = max(0.01, min(100.0, s))  # extreme 방지
        scale_x = scale_y = scale_z_val = s
    else:  # "none" 기본: 구조물 설계 크기 유지
        scale_x = scale_y = scale_z_val = 1.0

    # TIN 표면 Z (quad 중심에서 보간) — terrain 모드 공통
    z_surface = None
    if z_mode == "terrain":
        if pad_surface_z is not None:
            # 굴착 pad가 있으면 TIN 보간 대신 pad Z 직접 사용 (더 정확/안정)
            z_surface = float(pad_surface_z)
        elif terrain_mesh is not None:
            try:
                from scipy.interpolate import LinearNDInterpolator
                tin_pts = np.asarray(terrain_mesh.points, dtype=np.float64)
                interp = LinearNDInterpolator(tin_pts[:, :2], tin_pts[:, 2])
                zs = float(interp(q_center_local_x, q_center_local_y))
                if np.isnan(zs):
                    zs = float(np.median(tin_pts[:, 2]))
                z_surface = zs
            except Exception:
                z_surface = None

    out = []
    for mesh, color, opacity in filtered:
        try:
            pts = np.asarray(mesh.points, dtype=np.float64).copy()

            # 1) XY 중심을 원점으로 이동 (rotation/scale을 원점 기준으로)
            pts[:, 0] -= mesh_cx
            pts[:, 1] -= mesh_cy

            # 2) Anisotropic XY scale + Z scale (기하평균)
            pts[:, 0] *= scale_x
            pts[:, 1] *= scale_y
            if abs(scale_z_val - 1.0) > 1e-9:
                pts[:, 2] *= scale_z_val

            # 2.5) CW quad 보정: mesh Y 반전
            # 사용자가 시계방향으로 picks를 찍으면 quad의 local +Y가 edge 1→2 방향
            # (수학적 CCW 기준 -Y)이 됨. mesh의 +Y를 이쪽으로 정렬하려면 먼저 Y를
            # 반전한 후 표준 CCW 회전을 적용해야 함.
            if y_sign < 0:
                pts[:, 1] = -pts[:, 1]

            # 3) XY 회전 (상대 각도: TIN - detail - plan_frame_angle)
            x = pts[:, 0].copy()
            y = pts[:, 1].copy()
            pts[:, 0] = x * cos_a - y * sin_a
            pts[:, 1] = x * sin_a + y * cos_a

            # 4) XY를 quad 중심(TIN 로컬)으로 이동
            pts[:, 0] += q_center_local_x
            pts[:, 1] += q_center_local_y

            # 5) Z 이동 — 전체 구조물 최저 z_min을 TIN 표면에 맞춤
            # (컴포넌트 개별 z_min이 아니라 aggregate z_min 사용 → 상대 Z 보존)
            if z_mode == "terrain" and z_surface is not None:
                # embed_offset만큼 구조물을 아래로 내림 → TIN이 구조물 bottom보다
                # 살짝 위에 있어 아래에서 보면 TIN이 underside를 가림
                dz = z_surface - aggregate_z_min * scale_z_val - float(embed_offset)
                pts[:, 2] += dz
            elif z_mode == "offset":
                pts[:, 2] += float(z_offset)
            # z_mode == "design" 이면 Z 유지

            new_mesh = mesh.copy()
            new_mesh.points = pts
            out.append((new_mesh, color, opacity))
        except Exception:
            continue

    return out


def compute_orientation_from_points(points: list) -> float:
    """점들의 PCA로 주축 각도(도) 반환.

    Returns:
        주축과 X축 사이 각도 (-90 ~ +90 도)
    """
    if len(points) < 3:
        return 0.0
    arr = np.array(points, dtype=np.float64)
    arr = np.unique(arr, axis=0)
    if len(arr) < 3:
        return 0.0

    centered = arr[:, :2] - arr[:, :2].mean(axis=0)
    try:
        cov = np.cov(centered.T)
        eigenvals, eigenvecs = np.linalg.eigh(cov)
        # 최대 고유값 벡터가 주축
        idx = np.argmax(eigenvals)
        main_axis = eigenvecs[:, idx]
        angle = math.degrees(math.atan2(main_axis[1], main_axis[0]))
        # -90 ~ +90 범위로 정규화
        while angle > 90:
            angle -= 180
        while angle < -90:
            angle += 180
        return angle
    except Exception:
        return 0.0


def combine_meshes(mesh_groups: list[list[tuple]]) -> list[tuple]:
    """여러 구조물의 메쉬 리스트들을 하나로 합침."""
    combined = []
    for group in mesh_groups:
        combined.extend(group)
    return combined


if __name__ == "__main__":
    # 간단 테스트
    import pyvista as pv
    # 작은 박스
    box = pv.Box(bounds=(-1, 1, -1, 1, 0, 2))
    meshes = [(box, "#888888", 1.0)]
    placed = apply_placement(
        meshes=meshes,
        plan_centroid=(100.0, 50.0),
        rotation_deg=45.0,
        z_mode="offset",
        z_offset=10.0,
    )
    print(f"Original bounds: {meshes[0][0].bounds}")
    print(f"Placed bounds: {placed[0][0].bounds}")
    print("(X=99~101 → 99~101 회전, Y=49~51, Z=10~12)")