from __future__ import annotations import json from pathlib import Path from typing import Any import numpy as np from scipy import ndimage from scipy.interpolate import RectBivariateSpline, RBFInterpolator from skimage import measure # config 파라미터 로드 import config # 등고선 캐시 형식/추출 규칙이 바뀔 때 증가시킨다. CONTOUR_EXTRACTOR_VERSION = 3 def extract_contours_from_grid( x_coords: np.ndarray, y_coords: np.ndarray, z_grid: np.ndarray, valid_mask: np.ndarray, interval: float, min_interval: float = 0.5, scene_center: tuple[float, float, float] | None = None ) -> list[dict[str, Any]]: """ 정규 표고 격자로부터 등고선 라인을 추출합니다. scene_center가 (cx, cy, cz) 형태로 제공되면, 3D Scene 좌표계(Three.js)로 자동 변환하여 정점들을 저장합니다. """ interval = max(interval, min_interval) # z_grid의 유효값 범위 획득 finite_mask = np.isfinite(z_grid) if valid_mask is not None: finite_mask &= valid_mask if not finite_mask.any(): return [] z_min = float(np.min(z_grid[finite_mask])) z_max = float(np.max(z_grid[finite_mask])) # 등고선 생성 기준 레벨 결정 start_level = np.ceil(z_min / interval) * interval levels = np.arange(start_level, z_max, interval) if len(levels) == 0: return [] # 최대 레벨 제한 (과도한 루프 방지) if len(levels) > 500: new_interval = (z_max - z_min) / 100.0 levels = np.arange(np.ceil(z_min / new_interval) * new_interval, z_max, new_interval) interval = new_interval contours_geojson_list = [] # marching squares의 NaN 문제를 예방하기 위해, # 마스크가 아닌 영역(~valid_mask)을 z_min보다 충분히 낮은 sentinel(z_min - 1000)으로 채웁니다. z_grid_masked = z_grid.copy() if valid_mask is not None: z_grid_masked[~valid_mask] = z_min - 1000.0 # 유효한 영역 내의 NaN 값들도 함께 채워줍니다. invalid_mask = ~np.isfinite(z_grid_masked) if invalid_mask.any(): z_grid_masked[invalid_mask] = z_min - 1000.0 cx, cy, cz = scene_center if scene_center is not None else (0.0, 0.0, 0.0) # marching squares로 레벨별 등고선 추출 for level in levels: contours = measure.find_contours(z_grid_masked, level) for contour in contours: current_segment = [] for y_idx, x_idx in contour: # 격자 경계 클램프 x_idx_clamped = np.clip(x_idx, 0, len(x_coords) - 1) y_idx_clamped = np.clip(y_idx, 0, len(y_coords) - 1) # 인접 정수 인덱스 x_0 = int(np.floor(x_idx_clamped)) x_1 = int(np.ceil(x_idx_clamped)) y_0 = int(np.floor(y_idx_clamped)) y_1 = int(np.ceil(y_idx_clamped)) # 최종 정점이 유효 마스크 내부인지 엄격히 판별 is_valid = True if valid_mask is not None: if not (valid_mask[y_0, x_0] and valid_mask[y_0, x_1] and valid_mask[y_1, x_0] and valid_mask[y_1, x_1]): is_valid = False if not is_valid: # 유효하지 않은 영역을 만났다면, 현재까지 세그먼트가 존재할 경우 저장하고 새 세그먼트 준비 if len(current_segment) >= 2: mid_idx = len(current_segment) // 2 contours_geojson_list.append({ "level": float(level), "coordinates": current_segment, "label_position": current_segment[mid_idx] }) current_segment = [] continue # 소수부 비율 tx = x_idx_clamped - x_0 ty = y_idx_clamped - y_0 # 보간한 X, Y 좌표 계산 x_val = (1.0 - tx) * x_coords[x_0] + tx * x_coords[x_1] y_val = (1.0 - ty) * y_coords[y_0] + ty * y_coords[y_1] if scene_center is not None: # Three.js 씬 좌표 매핑 변환 적용 pt_x = round(float(x_val - cx), 3) pt_y = round(float(level - cz), 3) pt_z = round(float(-(y_val - cy)), 3) current_segment.append([pt_x, pt_y, pt_z]) else: current_segment.append([round(float(x_val), 3), round(float(y_val), 3), round(float(level), 3)]) # 마지막 남은 세그먼트가 유효하다면 추가 if len(current_segment) >= 2: mid_idx = len(current_segment) // 2 contours_geojson_list.append({ "level": float(level), "coordinates": current_segment, "label_position": current_segment[mid_idx] }) return contours_geojson_list def _load_footprint_mask( model_npz_path: Path, x_coords: np.ndarray, y_coords: np.ndarray, ) -> np.ndarray | None: """같은 source filter의 DTM valid_mask(실제 지면 footprint)를 현재 method 격자에 최근접 리샘플하여 반환한다. DTM 캐시가 없으면 None. (TIN/meshfree의 convex-hull 누출, NURBS/implicit의 직사각형 영역 누출을 모두 차단) """ stem = Path(model_npz_path).stem # 예: "tin_grid_min_z_smooth" 또는 "tin_grid_min_z" if stem.endswith("_smooth"): stem = stem[:-7] parts = stem.split("_", 1) if len(parts) < 2: return None filter_key = parts[1] # "grid_min_z", "csf", "pmf" dtm_path = Path(model_npz_path).parent / f"dtm_{filter_key}.npz" if not dtm_path.exists(): return None try: d = np.load(dtm_path) dtm_x = np.asarray(d["x"]).ravel() dtm_y = np.asarray(d["y"]).ravel() dtm_mask = np.asarray(d["valid_mask"], dtype=bool) except Exception: return None if len(dtm_x) < 2 or len(dtm_y) < 2: return None def _nearest_idx(axis: np.ndarray, coords: np.ndarray) -> np.ndarray: ascending = bool(axis[0] <= axis[-1]) a = axis if ascending else axis[::-1] idx = np.clip(np.searchsorted(a, coords), 1, len(a) - 1) idx = np.where(np.abs(a[idx - 1] - coords) <= np.abs(a[idx] - coords), idx - 1, idx) return idx if ascending else (len(axis) - 1 - idx) xi = _nearest_idx(dtm_x, np.asarray(x_coords, dtype=np.float64)) yi = _nearest_idx(dtm_y, np.asarray(y_coords, dtype=np.float64)) return dtm_mask[np.ix_(yi, xi)] def _apply_footprint( model_npz_path: Path, x_coords: np.ndarray, y_coords: np.ndarray, valid_mask: np.ndarray, ) -> np.ndarray: """현재 valid_mask에 DTM 기반 footprint 마스크를 교집합으로 적용한다. 형상이 약간 다르면 zoom/nearest로 자동 형상 조절하여 강제 적용한다.""" fp = _load_footprint_mask(model_npz_path, x_coords, y_coords) if fp is not None: if fp.shape == valid_mask.shape: return valid_mask & fp else: # 해상도가 다를 경우 최근접 리샘플링하여 크기를 맞춤 from scipy.ndimage import zoom zoom_y = valid_mask.shape[0] / fp.shape[0] zoom_x = valid_mask.shape[1] / fp.shape[1] fp_resized = zoom(fp.astype(float), (zoom_y, zoom_x), order=0) > 0.5 if fp_resized.shape == valid_mask.shape: return valid_mask & fp_resized return valid_mask def _tin_face_coverage_mask( vertices: np.ndarray, faces: np.ndarray, xx: np.ndarray, yy: np.ndarray, ) -> np.ndarray: """저장된 TIN 면이 실제로 덮는 XY 영역만 True로 반환한다. 정점으로 Delaunay를 다시 수행하면 TIN 생성 단계에서 제거한 긴 면과 footprint 경계 면이 되살아날 수 있다. 따라서 실제 저장된 faces의 경계 에지를 polygonize하여 렌더링 메시와 동일한 coverage를 사용한다. """ vertices = np.asarray(vertices) faces = np.asarray(faces, dtype=np.int64) if vertices.ndim != 2 or vertices.shape[1] < 2 or not len(faces): return np.zeros(xx.shape, dtype=bool) # 삼각형마다 세 변을 모은 뒤 한 번만 나타나는 변을 메시 경계로 본다. edges = np.vstack((faces[:, [0, 1]], faces[:, [1, 2]], faces[:, [2, 0]])) edges = np.sort(edges, axis=1) unique_edges, counts = np.unique(edges, axis=0, return_counts=True) boundary_edges = unique_edges[counts == 1] if not len(boundary_edges): return np.zeros(xx.shape, dtype=bool) from shapely import get_parts, intersects_xy, linestrings, polygonize, union_all boundary_lines = linestrings(vertices[boundary_edges, :2]) polygons = list(get_parts(polygonize(boundary_lines))) if not polygons: return np.zeros(xx.shape, dtype=bool) coverage = union_all(polygons) # shapely.intersects_xy의 x, y 인자는 1차원 float64 배열이어야 하므로 평탄화 후 타입 캐스팅 적용 xx_flat = np.asarray(xx, dtype=np.float64).ravel() yy_flat = np.asarray(yy, dtype=np.float64).ravel() res_flat = np.asarray(intersects_xy(coverage, xx_flat, yy_flat), dtype=bool) return res_flat.reshape(xx.shape) def extract_contours( model_npz_path: Path, representation: str, interval: float, target_grid_m: float = 1.0, scene_center: tuple[float, float, float] | None = None ) -> list[dict[str, Any]]: """ 15종 표현형식의 npz 모델 파일로부터 표고 격자를 환원한 후, 정해진 interval 간격의 등고선 리스트를 추출합니다. """ model_npz_path = Path(model_npz_path) if not model_npz_path.exists(): raise FileNotFoundError(f"모델 파일이 존재하지 않습니다: {model_npz_path}") data = np.load(model_npz_path) # 1. 표현형식별 정규 격자(elevation grid) 환원 if representation == "regular_grid": # DTM 격자 다운샘플 x_coords = data["x"] y_coords = data["y"] z_grid = data["z"] valid_mask = data["valid_mask"] # dtm의 dtm_grid_resolution_meters는 보통 0.05m이므로 1m 등으로 다운샘플 current_res = (x_coords[-1] - x_coords[0]) / (len(x_coords) - 1) step = max(1, int(round(target_grid_m / current_res))) if step > 1: x_coords_sub = x_coords[::step] y_coords_sub = y_coords[::step] z_grid_sub = z_grid[::step, ::step] valid_mask_sub = valid_mask[::step, ::step] return extract_contours_from_grid(x_coords_sub, y_coords_sub, z_grid_sub, valid_mask_sub, interval, scene_center=scene_center) else: return extract_contours_from_grid(x_coords, y_coords, z_grid, valid_mask, interval, scene_center=scene_center) elif representation == "triangular_mesh": # TIN 메쉬 보간 vertices = data["vertices"] # (N, 3) faces = data["faces"] # 1.0m 격자로 샘플링 x_min, x_max = float(np.min(vertices[:, 0])), float(np.max(vertices[:, 0])) y_min, y_max = float(np.min(vertices[:, 1])), float(np.max(vertices[:, 1])) cols = max(2, int(np.ceil((x_max - x_min) / target_grid_m)) + 1) rows = max(2, int(np.ceil((y_max - y_min) / target_grid_m)) + 1) x_coords = np.linspace(x_min, x_max, cols, dtype=np.float32) y_coords = np.linspace(y_min, y_max, rows, dtype=np.float32) xx, yy = np.meshgrid(x_coords, y_coords) # scipy.interpolate.griddata를 사용해 TIN 지면 보간 from scipy.interpolate import griddata z_grid = griddata(vertices[:, :2], vertices[:, 2], (xx, yy), method="linear") # 중요: vertices를 다시 Delaunay하면 생성 단계에서 제거한 삼각형이 # 복원된다. 저장된 faces의 실제 XY coverage만 유효 영역으로 사용한다. face_mask = _tin_face_coverage_mask(vertices, faces, xx, yy) valid_mask = np.isfinite(z_grid) & face_mask # 실제 지면 footprint와 교집합 (DTM 기준) valid_mask = _apply_footprint(model_npz_path, x_coords, y_coords, valid_mask) return extract_contours_from_grid(x_coords, y_coords, z_grid, valid_mask, interval, scene_center=scene_center) elif representation == "bspline_surface": # NURBS 제어점 로드 및 평가 control_x = data["control_x"] control_y = data["control_y"] control_z = data["control_z"] degree = int(data["degree"][0]) spline = RectBivariateSpline( control_y, control_x, control_z, kx=min(degree, len(control_y) - 1), ky=min(degree, len(control_x) - 1), s=float(len(control_x) * len(control_y)) * 0.01 ) # 1.0m 격자로 평가 x_min, x_max = float(control_x[0]), float(control_x[-1]) y_min, y_max = float(control_y[0]), float(control_y[-1]) cols = max(2, int(np.ceil((x_max - x_min) / target_grid_m)) + 1) rows = max(2, int(np.ceil((y_max - y_min) / target_grid_m)) + 1) x_coords = np.linspace(x_min, x_max, cols, dtype=np.float32) y_coords = np.linspace(y_min, y_max, rows, dtype=np.float32) z_grid = np.asarray(spline(y_coords, x_coords), dtype=np.float32) # NURBS는 직사각형 도메인 전체를 평가하므로 실제 지면 footprint로 제한한다. valid_mask = np.ones_like(z_grid, dtype=bool) valid_mask = _apply_footprint(model_npz_path, x_coords, y_coords, valid_mask) return extract_contours_from_grid(x_coords, y_coords, z_grid, valid_mask, interval, scene_center=scene_center) elif representation == "local_rbf_height_field": # Implicit RBF 평가 centers_xy = data["centers_xy"] center_z = data["center_z"] smoothing = float(data["smoothing"][0]) interpolator = RBFInterpolator( centers_xy.astype(np.float64), center_z.astype(np.float64), neighbors=min(64, len(centers_xy)), smoothing=smoothing, kernel="thin_plate_spline" ) x_min, x_max = float(np.min(centers_xy[:, 0])), float(np.max(centers_xy[:, 0])) y_min, y_max = float(np.min(centers_xy[:, 1])), float(np.max(centers_xy[:, 1])) cols = max(2, int(np.ceil((x_max - x_min) / target_grid_m)) + 1) rows = max(2, int(np.ceil((y_max - y_min) / target_grid_m)) + 1) x_coords = np.linspace(x_min, x_max, cols, dtype=np.float32) y_coords = np.linspace(y_min, y_max, rows, dtype=np.float32) xx, yy = np.meshgrid(x_coords, y_coords) query = np.column_stack([xx.ravel(), yy.ravel()]) z_values = interpolator(query).astype(np.float32) z_grid = z_values.reshape(len(y_coords), len(x_coords)) # Implicit RBF도 직사각형 도메인 전체를 평가하므로 실제 지면 footprint로 제한한다. valid_mask = np.ones_like(z_grid, dtype=bool) valid_mask = _apply_footprint(model_npz_path, x_coords, y_coords, valid_mask) return extract_contours_from_grid(x_coords, y_coords, z_grid, valid_mask, interval, scene_center=scene_center) elif representation == "meshfree_surfels": # Meshfree 포인트 데이터 보간 points = data["points"] # (N, 3) x_min, x_max = float(np.min(points[:, 0])), float(np.max(points[:, 0])) y_min, y_max = float(np.min(points[:, 1])), float(np.max(points[:, 1])) cols = max(2, int(np.ceil((x_max - x_min) / target_grid_m)) + 1) rows = max(2, int(np.ceil((y_max - y_min) / target_grid_m)) + 1) x_coords = np.linspace(x_min, x_max, cols, dtype=np.float32) y_coords = np.linspace(y_min, y_max, rows, dtype=np.float32) xx, yy = np.meshgrid(x_coords, y_coords) from scipy.interpolate import griddata from scipy.spatial import Delaunay z_grid = griddata(points[:, :2], points[:, 2], (xx, yy), method="linear") valid_mask = np.isfinite(z_grid) # Meshfree도 griddata(linear)라 TIN과 동일하게 convex hull 경계 마스킹 try: tri = Delaunay(points[:, :2]) query_pts = np.column_stack([xx.ravel(), yy.ravel()]) hull_inside = tri.find_simplex(query_pts) >= 0 hull_mask = hull_inside.reshape(xx.shape) valid_mask = valid_mask & hull_mask except Exception: pass valid_mask = _apply_footprint(model_npz_path, x_coords, y_coords, valid_mask) return extract_contours_from_grid(x_coords, y_coords, z_grid, valid_mask, interval, scene_center=scene_center) else: raise ValueError(f"지원하지 않는 표현 방식입니다: {representation}")