diff --git a/B04_PreProcess/B04_PreProcess_Engine_Contour.py b/B04_PreProcess/B04_PreProcess_Engine_Contour.py index 3ed70a52..7ede4004 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_Contour.py +++ b/B04_PreProcess/B04_PreProcess_Engine_Contour.py @@ -10,7 +10,7 @@ from pathlib import Path from typing import Any import numpy as np -from scipy.interpolate import RBFInterpolator, RectBivariateSpline +from scipy.interpolate import RBFInterpolator from skimage import measure # 등고선 캐시 형식/추출 규칙이 바뀔 때 증가시킨다. @@ -36,14 +36,17 @@ def extract_contours_from_grid( z_min = float(np.min(z_grid[finite_mask])) z_max = float(np.max(z_grid[finite_mask])) + span = z_max - z_min + if span <= 0: + return [] + # 레벨 수 상한은 arange **앞에서** 건다. 뒤에 두면 표고 범위가 비정상적으로 넓을 때 + # np.arange가 배열을 만들다 "Maximum allowed size exceeded"로 먼저 터진다. + if span / interval > 500: + interval = span / 100.0 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: list[dict[str, Any]] = [] @@ -398,16 +401,15 @@ def extract_contours( ) if representation == "bspline_surface": + # 모델 빌더와 같은 곡면을 다시 만든다 — 평활 계수·발산 방어를 공유한다. + from B04_PreProcess.B04_PreProcess_Engine_ModelContext import ( + evaluate_nurbs_spline, + fit_nurbs_spline, + ) + control_x, control_y, control_z = data["control_x"], data["control_y"], 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, - ) + spline, z_range = fit_nurbs_spline(control_x, control_y, control_z, degree) x_coords, y_coords = _grid_axes( float(control_x[0]), float(control_x[-1]), @@ -415,7 +417,9 @@ def extract_contours( float(control_y[-1]), target_grid_m, ) - z_grid = np.asarray(spline(y_coords, x_coords), dtype=np.float32) + z_grid = evaluate_nurbs_spline( + spline, y_coords, x_coords, z_range, Path(model_npz_path).stem + ).astype(np.float32) valid_mask = _apply_footprint( model_npz_path, x_coords, y_coords, np.ones_like(z_grid, dtype=bool) ) diff --git a/B04_PreProcess/B04_PreProcess_Engine_ModelBuild.py b/B04_PreProcess/B04_PreProcess_Engine_ModelBuild.py index 09f032c8..baa33c69 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_ModelBuild.py +++ b/B04_PreProcess/B04_PreProcess_Engine_ModelBuild.py @@ -8,7 +8,7 @@ from pathlib import Path from typing import Any import numpy as np -from scipy.interpolate import RBFInterpolator, RectBivariateSpline +from scipy.interpolate import RBFInterpolator from scipy.spatial import Delaunay from B04_PreProcess.B04_PreProcess_Engine_ModelContext import ( @@ -17,6 +17,8 @@ from B04_PreProcess.B04_PreProcess_Engine_ModelContext import ( artifact_size, atomic_npz, clip_and_compact_mesh, + evaluate_nurbs_spline, + fit_nurbs_spline, grid_faces, grid_vertices, with_footprint, @@ -125,18 +127,13 @@ def build_nurbs( control_resolution = max(patch_size / max(controls - 1, 1), 0.25) x_control, y_control, z_control = context.grid(control_resolution) progress(30) - spline = RectBivariateSpline( - y_control, - x_control, - z_control, - kx=min(degree, len(y_control) - 1), - ky=min(degree, len(x_control) - 1), - s=float(len(x_control) * len(y_control)) * 0.01, - ) + spline, z_range = fit_nurbs_spline(x_control, y_control, z_control, degree) x_preview, y_preview, _ = context.preview_grid( float(context.config["dtm_grid_resolution_meters"]) ) - z_preview = np.asarray(spline(y_preview, x_preview), dtype=np.float32) + z_preview = evaluate_nurbs_spline(spline, y_preview, x_preview, z_range, stem).astype( + np.float32 + ) progress(65) vertices = grid_vertices(x_preview, y_preview, z_preview) faces = grid_faces(len(y_preview), len(x_preview)) diff --git a/B04_PreProcess/B04_PreProcess_Engine_ModelContext.py b/B04_PreProcess/B04_PreProcess_Engine_ModelContext.py index a113f652..e58fe9e4 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_ModelContext.py +++ b/B04_PreProcess/B04_PreProcess_Engine_ModelContext.py @@ -7,6 +7,7 @@ GLB/PLY 프리뷰 및 npz 모델을 원자적으로 저장하는 공통 기능 import hashlib import json +import logging import math from dataclasses import dataclass, field from pathlib import Path @@ -15,9 +16,12 @@ from typing import Any, Callable import numpy as np import trimesh from scipy import ndimage +from scipy.interpolate import RectBivariateSpline from common_util.common_util_atomic import atomic_write_bytes, atomic_write_npz +logger = logging.getLogger(__name__) + MODEL_VERSION = 1 MODEL_METHODS = ("tin", "dtm", "nurbs", "implicit", "meshfree") SOURCE_FILTER_LABELS = {"grid_min_z": "Grid Min-Z", "csf": "CSF", "pmf": "PMF"} @@ -26,6 +30,62 @@ ProgressCallback = Callable[[int], None] # 대용량 포인트 배치 처리 크기 _BATCH_SIZE = 500_000 +# NURBS 평활 계수 — 제어점당 허용 잔차제곱(m²). +# 0.01(RMS 0.1m)은 요구가 너무 빡빡해 FITPACK이 제어점 수(122x113)보다 많은 +# knot(126x117)을 밀어넣고 계수가 발산했다 — csf 표고 1e105, grid_min_z 1e10 +# (2026-08-17 실측). 1.0(RMS 1m)이면 knot이 17~35로 떨어지고 세 지면 필터 모두 +# 데이터 표고 범위 안에 머문다. NURBS는 평활 곡면 표현이라 이 정도 완화가 맞다. +NURBS_RESIDUAL_PER_CONTROL_POINT = 1.0 + + +def fit_nurbs_spline( + x_control: np.ndarray, y_control: np.ndarray, z_control: np.ndarray, degree: int +) -> tuple[RectBivariateSpline, tuple[float, float]]: + """제어 격자에 B-spline 곡면을 맞추고, 발산 판정용 표고 허용 범위를 함께 준다. + + 모델 빌더와 등고선 엔진이 같은 곡면을 각각 만들므로 여기서 한 번에 정의한다. + """ + spline = RectBivariateSpline( + y_control, + x_control, + z_control, + kx=min(degree, len(y_control) - 1), + ky=min(degree, len(x_control) - 1), + s=float(len(x_control) * len(y_control)) * NURBS_RESIDUAL_PER_CONTROL_POINT, + ) + z_min = float(np.min(z_control)) + z_max = float(np.max(z_control)) + margin = max((z_max - z_min) * 0.5, 1.0) + return spline, (z_min - margin, z_max + margin) + + +def evaluate_nurbs_spline( + spline: RectBivariateSpline, + y_coords: np.ndarray, + x_coords: np.ndarray, + z_range: tuple[float, float], + label: str = "", +) -> np.ndarray: + """스플라인을 평가하고 표고 허용 범위를 벗어난 값을 잘라낸다. + + 평활 스플라인은 제어점이 성긴 구석에서 발산할 수 있다. 그대로 두면 float32 + 캐스팅에서 inf가 되고 프리뷰 색상(nan)과 등고선 레벨 산출(`np.arange`가 + "Maximum allowed size exceeded"로 실패)까지 번진다. + """ + z_values = np.asarray(spline(y_coords, x_coords), dtype=np.float64) + low, high = z_range + outside = np.count_nonzero(~np.isfinite(z_values) | (z_values < low) | (z_values > high)) + if outside: + logger.warning( + "NURBS 곡면이 표고 범위를 벗어나 잘라냈습니다: %s%d개 셀 (허용 %.1f~%.1fm)", + f"{label} " if label else "", + int(outside), + low, + high, + ) + z_values = np.nan_to_num(z_values, nan=low, posinf=high, neginf=low) + return np.clip(z_values, low, high) + def config_signature(config: dict[str, Any]) -> str: """지오메트리 원본과 무관한 등고선·스무딩 설정을 제외한 캐시 서명."""