diff --git a/B04_PreProcess/B04_PreProcess_Engine_SheetSurface.py b/B04_PreProcess/B04_PreProcess_Engine_SheetSurface.py index a5f57920..e475fbc3 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_SheetSurface.py +++ b/B04_PreProcess/B04_PreProcess_Engine_SheetSurface.py @@ -34,6 +34,7 @@ from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import ( from config.config_system import ( SHEET_SURFACE_GRID_M, SHEET_SURFACE_MARGIN_M, + SHEET_SURFACE_MIDLINE_MAX_M, SURFACE_MAX_PREVIEW_VERTICES, ) @@ -101,6 +102,68 @@ def _preview_mesh( return clip_and_compact_mesh(vertices, faces, pv.reshape(-1)) +def _densify_between_contours(spec: Any, features: list[dict[str, Any]], cloud: Any) -> Any: + """인접 표고 등고선 쌍의 등거리 중간선을 중간 표고 정점으로 추가한다. + + 등고선 정점만으로 Delaunay TIN을 치면 같은 표고 정점 3개짜리 평탄 삼각형이 + 굴곡부·능선에서 계단(terrace)을 만든다. 두 등고선의 "사이점"(등거리선)에 + 평균 표고 정점을 넣어 삼각형이 반드시 서로 다른 표고를 잇게 한다 + (2026-08-30 사용자 지적). 최고 등고선 안쪽 봉우리·같은 표고 사이 골짜기 + 바닥은 짝이 없어 평탄으로 남는다 — 원천 데이터 한계. + """ + from scipy.ndimage import distance_transform_edt + + from B04_PreProcess.B04_PreProcess_Engine_Watershed_Descent import rasterize_contours + from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import ContourCloud + + burned, levels = rasterize_contours(spec, features, None) + present = [level for level in sorted(levels) if bool((burned == level).any())] + if len(present) < 2: + return cloud + + def _distance(level: float) -> np.ndarray: + return distance_transform_edt(burned != level, sampling=spec.cell_m).astype(np.float32) + + # 1차: 셀별 최근접 등고선 레벨 — 등거리 조건은 멀리 있는 다른 표고 쌍에도 우연히 + # 성립하므로(예: 525~530 골짜기에 555/560 중간점이 박힘), "그 지점의 최근접 + # 등고선이 바로 그 쌍"일 때만 사이점으로 인정한다. + best_distance = np.full(burned.shape, np.inf, dtype=np.float32) + best_index = np.full(burned.shape, -1, dtype=np.int16) + for index, level in enumerate(present): + distance = _distance(level) + take = distance < best_distance + best_distance[take] = distance[take] + best_index[take] = index + + xs = spec.cell_centers_x() + ys = spec.cell_centers_y() + extra_xy: list[np.ndarray] = [] + extra_z: list[np.ndarray] = [] + previous_distance = _distance(present[0]) + for index in range(1, len(present)): + distance = _distance(present[index]) + near = (previous_distance < SHEET_SURFACE_MIDLINE_MAX_M) & ( + distance < SHEET_SURFACE_MIDLINE_MAX_M + ) + pair_is_nearest = (best_index == index - 1) | (best_index == index) + midline = near & pair_is_nearest & (np.abs(previous_distance - distance) <= spec.cell_m) + rows, cols = np.nonzero(midline) + if len(rows): + # 등고선 정점 재샘플 밀도에 맞춰 성기게 담는다 — TIN 비용 절제. + rows, cols = rows[::5], cols[::5] + extra_xy.append(np.column_stack([xs[cols], ys[rows]])) + extra_z.append(np.full(len(rows), (present[index - 1] + present[index]) / 2.0)) + previous_distance = distance + if not extra_xy: + return cloud + added = int(sum(len(a) for a in extra_z)) + logger.info("도엽 서피스: 중간 보간선 정점 %d개 추가 (등고 %d단)", added, len(present)) + return ContourCloud( + xy=np.vstack([cloud.xy, *extra_xy]), + z=np.concatenate([cloud.z, *extra_z]), + ) + + def build_sheet_surface_model( project_root: Path, processed_dir: Path, @@ -130,6 +193,8 @@ def build_sheet_surface_model( return None spec = grid_spec_from_bounds(x_min, y_min, x_max, y_max, SHEET_SURFACE_GRID_M) + # 계단(평탄 삼각형) 해소 — 인접 등고선 사이 중간 보간선을 정점으로 추가. + cloud = _densify_between_contours(spec, features, cloud) surface = interpolate_elevation(spec, cloud) # (R, C), 북→남 행 순서, 외부 NaN # DtmGridSampler 규약에 맞춰 y 오름차순으로 뒤집어 저장한다. diff --git a/config/config_system.py b/config/config_system.py index e49ca9e6..8272c5ac 100644 --- a/config/config_system.py +++ b/config/config_system.py @@ -188,6 +188,9 @@ SURFACE_CONTOUR_GRID_RESOLUTION_M = float(os.getenv("SURFACE_CONTOUR_GRID_RESOLU SHEET_SURFACE_MARGIN_M = float(os.getenv("SHEET_SURFACE_MARGIN_M", "300.0")) # 도엽등고선 DTM 격자 한 변(m). LAS DTM·등고선 캐시와 같은 1m(사용자 확정값). SHEET_SURFACE_GRID_M = float(os.getenv("SHEET_SURFACE_GRID_M", "1.0")) +# 중간 보간선(인접 등고선 쌍의 등거리선) 탐색 최대 거리(m). 등고선 간격이 이보다 넓은 +# 완경사면은 보간선을 넣지 않는다. 5m 주곡선 기준 경사 2.5%까지 커버. +SHEET_SURFACE_MIDLINE_MAX_M = float(os.getenv("SHEET_SURFACE_MIDLINE_MAX_M", "200.0")) # 일반 사용자 WF1 자동 확정 기본값 SURFACE_CONFIRM_DEFAULT_FILTER = os.getenv("SURFACE_CONFIRM_DEFAULT_FILTER", "csf")