diff --git a/B04_PreProcess/B04_PreProcess_Engine_SheetSurface.py b/B04_PreProcess/B04_PreProcess_Engine_SheetSurface.py index fd8932f4..9b215006 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_SheetSurface.py +++ b/B04_PreProcess/B04_PreProcess_Engine_SheetSurface.py @@ -9,8 +9,9 @@ LAS 없는 설계(2026-08-30 사용자 확정)의 지형 원천이자, LAS가 격자는 `SHEET_SURFACE_GRID_M`(1m). 순서는 **2D 먼저, 메시는 맨 마지막**이다(2026-08-30 사용자 지시): - ① 등고 라인을 격자에 굽고 ② 계곡 구조선 앵커를 제약으로 더한 뒤 - ③ 등고선 사이 거리 비례 보간으로 표고 격자를 만들고 ④ 마루를 캡한다. + ① 등고 라인을 격자에 굽고 ② 등고선 사이 거리 비례 보간으로 표고 격자를 만든 뒤 + ③ 폐합 등고선 안쪽(마루·웅덩이)을 바깥 사면 경사로 연장하고 ④ 라인을 고정한 채 + 완화(라플라스)해 등고 간격을 고르게 한다. 격자에서 뽑는 1m 등고선이 곧 2D 보간선이며, 메시(glb)는 그 격자의 표현일 뿐이다. """ @@ -34,6 +35,7 @@ from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import grid_spec_from_b from config.config_system import ( SHEET_SURFACE_GRID_M, SHEET_SURFACE_MARGIN_M, + SHEET_SURFACE_RELAX_ITERATIONS, SURFACE_MAX_PREVIEW_VERTICES, ) @@ -41,7 +43,6 @@ logger = logging.getLogger(__name__) # 도엽 병합 산출물 파일명 (B04_PreProcess_Router_Watershed와 같은 값) _CONTOUR_FILE = "도엽_등고선.geojson" -_STREAM_FILE = "도엽_하천중심선.geojson" # 산출 모델 식별자 — surface_models.generation_params.source_filter 및 파일 stem에 쓴다. SHEET_SOURCE_FILTER = "sheet" @@ -104,72 +105,107 @@ def _preview_mesh( return clip_and_compact_mesh(vertices, faces, pv.reshape(-1)) -def _stream_breakline_vertices( - spec: Any, burned: np.ndarray, stream_features: list[dict[str, Any]] -) -> tuple[np.ndarray, np.ndarray] | None: - """계곡 기준선(하천중심선)을 따라 등고 교차점 사이를 보간한 정점열을 만든다. +def _rasterize_contour_levels(spec: Any, features: list[dict[str, Any]]) -> np.ndarray: + """등고 라인을 격자에 굽는다 — 셀 = 그 위를 지나는 라인의 표고, 그 외 NaN. - 구조선 기반 보간(2026-08-30 사용자 확정): 계곡 축을 먼저 세우고, 축이 등고선과 - 만나는 점을 Z 앵커로 삼아 앵커 사이를 선 길이 비례로 보간한다. 계곡 바닥이 - 등고선 사이에서도 연속으로 내려가는 가상 종단이 되어 골짜기 평탄·역경사가 준다. + 배수유역의 `rasterize_contours()`와 두 가지가 다르다(둘 다 서피스 품질 때문이다): + · `all_touched=False` — 스치는 셀까지 칠하면 라인이 2px 두께가 되고, 그 폭만큼 + 정확히 등고 표고인 **평탄 띠**가 생겨 사이 1m 등고선 간격이 찌그러진다 + (2026-08-30 사용자 지적: 보간선이 등간격이 아님). + · 길이 필터 없음 — 봉우리 폐합 링 같은 짧은 등고선을 버리면 그 일대가 통째로 + 평평해진다. 배수유역은 노이즈를 버려야 하지만 지형면은 다 있어야 한다. + + 같은 셀을 두 표고가 지나면 낮은 쪽을 남긴다(배수유역과 같은 규칙). """ - from shapely import segmentize + from rasterio.features import rasterize from shapely.geometry import shape - from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import iter_linestrings + from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import ( + ELEVATION_KEYS, + grid_transform, + iter_linestrings, + ) - sample_step_m = 2.0 - vertex_stride = 3 # 2m 샘플 → 6m 간격 정점 (등고 정점 재샘플 밀도와 유사) - collected_xy: list[np.ndarray] = [] - collected_z: list[np.ndarray] = [] - for feature in stream_features: + by_level: dict[float, list[Any]] = {} + for feature in features: geometry = feature.get("geometry") if not geometry: continue + properties = feature.get("properties") or {} + elevation = next( + (float(properties[key]) for key in ELEVATION_KEYS if properties.get(key) is not None), + None, + ) + if elevation is None: + continue try: parsed = shape(geometry) except Exception: # noqa: BLE001 continue for line in iter_linestrings(parsed): - coords = np.asarray(segmentize(line, sample_step_m).coords, dtype=np.float64) - if len(coords) < 3: - continue - xy = coords[:, :2] - arc = np.concatenate([[0.0], np.cumsum(np.hypot(*np.diff(xy, axis=0).T))]) - row, col = spec.world_to_rc(xy[:, 0].copy(), xy[:, 1].copy()) - inside = (row >= 0) & (col >= 0) - levels_at = np.full(len(xy), np.nan, dtype=np.float32) - levels_at[inside] = burned[row[inside], col[inside]] - hit = np.isfinite(levels_at) - if hit.sum() < 2: - continue - # 등고선 통과 구간(연속 같은 표고 샘플 묶음) 하나 = 앵커 하나. - anchor_arc: list[float] = [] - anchor_z: list[float] = [] - run_start: int | None = None - for i in range(len(xy) + 1): - is_hit = i < len(xy) and hit[i] - if is_hit and run_start is not None and levels_at[i] != levels_at[run_start]: - is_hit = False # 표고가 바뀌면 묶음을 끊는다 - if is_hit: - if run_start is None: - run_start = i - elif run_start is not None: - anchor_arc.append(float(arc[run_start : i if i <= len(xy) else len(xy)].mean())) - anchor_z.append(float(levels_at[run_start])) - run_start = i if i < len(xy) and hit[i] else None - if len(anchor_arc) < 2: - continue - # 앵커 사이 구간만 채택 — 첫 앵커 이전·마지막 앵커 이후는 근거가 없다. - span = (arc >= anchor_arc[0]) & (arc <= anchor_arc[-1]) & inside & ~hit - take = np.flatnonzero(span)[::vertex_stride] - if not len(take): - continue - collected_xy.append(xy[take]) - collected_z.append(np.interp(arc[take], anchor_arc, anchor_z)) - if not collected_xy: - return None - return np.vstack(collected_xy), np.concatenate(collected_z) + by_level.setdefault(elevation, []).append(line) + + burned = np.full((spec.n_rows, spec.n_cols), np.nan, dtype=np.float32) + transform = grid_transform(spec) + for elevation in sorted(by_level, reverse=True): + stamp = rasterize( + [(line, 1) for line in by_level[elevation]], + out_shape=(spec.n_rows, spec.n_cols), + transform=transform, + fill=0, + dtype="uint8", + all_touched=False, + ).astype(bool) + burned[stamp] = elevation + logger.info( + "도엽 서피스: 등고 라인 %d단을 격자에 굽어 %d셀", + len(by_level), + int(np.isfinite(burned).sum()), + ) + return burned + + +def _relax_surface(surface: np.ndarray, fixed: np.ndarray, iterations: int) -> None: + """등고 라인을 고정한 채 이웃 평균으로 표고를 **살짝만** 다듬는다 (in-place). + + 거리 보간은 두 라인의 굽은 정도가 다르면 능선·계곡 중심축에 각진 자국이 남는다. + 이웃 평균 몇 회로 그 자국만 지운다. + + **끝까지 수렴시키면 안 된다.** 수렴한 해는 라플라스 방정식의 해(harmonic)인데, + 원뿔·능선 형상 z=r은 harmonic이 아니라 biharmonic이라 harmonic 해로 끌고 가면 + 마루가 눌리고 평탄 셀이 늘어난다. 등고선 DEM의 표준인 ANUDEM/Topo to Raster가 + 라플라스가 아니라 **thin plate spline(biharmonic)** 을 쓰는 이유가 이것이다 + (Hutchinson 1988/89). 거리 보간은 원뿔을 정확히 재현하므로 그것을 바탕으로 두고 + 다듬기만 한다. + + 합성 원뿔 검증(2026-08-30): 거리보간만 오차 0.121m·간격 CV 2.18 → + 5회 다듬기 0.118m·1.43(개선) → 20회 이상 CV 50↑·평탄 셀 증가(악화). + """ + if iterations <= 0: + return + free = ~fixed & np.isfinite(surface) + if not free.any(): + return + # 체커보드(red-black) 순서 — 한 색을 갱신할 때 이웃(다른 색)은 그대로라 제자리 + # 갱신이 안전하다. 과완화는 쓰지 않는다(omega=1): 빨리 수렴시킬수록 harmonic 해에 + # 가까워져 능선이 눌린다. 여기서 원하는 건 수렴이 아니라 자국 제거다. + omega = np.float32(1.0) + rows, cols = np.indices(surface.shape) + red = free & (((rows + cols) & 1) == 0) + black = free & ~red + padded = np.zeros((surface.shape[0] + 2, surface.shape[1] + 2), dtype=np.float32) + for _ in range(iterations): + for colour in (red, black): + padded[1:-1, 1:-1] = surface + padded[0, 1:-1] = surface[0] + padded[-1, 1:-1] = surface[-1] + padded[1:-1, 0] = surface[:, 0] + padded[1:-1, -1] = surface[:, -1] + neighbours = ( + padded[:-2, 1:-1] + padded[2:, 1:-1] + padded[1:-1, :-2] + padded[1:-1, 2:] + ) * np.float32(0.25) + surface[colour] += omega * (neighbours[colour] - surface[colour]) + logger.info("도엽 서피스: 등고선 고정 완화 %d회 (자유 셀 %d개)", iterations, int(free.sum())) def _interpolate_between_contours(burned: np.ndarray, cell_m: float) -> np.ndarray: @@ -229,27 +265,9 @@ def _interpolate_between_contours(burned: np.ndarray, cell_m: float) -> np.ndarr return surface -def _burn_stream_anchors( - spec: Any, burned: np.ndarray, stream_features: list[dict[str, Any]] -) -> int: - """계곡 구조선 앵커 보간값을 등고 격자에 제약으로 굽는다. 구운 셀 수 반환.""" - vertices = _stream_breakline_vertices(spec, burned, stream_features) - if vertices is None: - return 0 - xy, z = vertices - row, col = spec.world_to_rc(xy[:, 0].copy(), xy[:, 1].copy()) - inside = (row >= 0) & (col >= 0) - row, col, z = row[inside], col[inside], z[inside] - free = ~np.isfinite(burned[row, col]) # 등고 라인 위는 덮지 않는다 - # 1m로 양자화 — 거리장을 표고별로 한 번씩 도는 방식이라 레벨 수가 곧 비용이다. - # 원천이 5m 주곡선이므로 0.5m 이내 반올림은 정확도에 영향이 없다. - burned[row[free], col[free]] = np.round(z[free]).astype(np.float32) - return int(free.sum()) - - def _resolve_enclosed_interiors( burned: np.ndarray, present: list[float], surface: np.ndarray, cell_m: float, interval_m: float -) -> int: +) -> np.ndarray: """폐합 등고선 안쪽(봉우리·웅덩이)을 바깥 사면 경사로 연장한다 (in-place). 거리 보간은 "가장 가까운 서로 다른 두 라인 사이"를 채우므로, 마지막 등고선 @@ -260,7 +278,7 @@ def _resolve_enclosed_interiors( 폐합 라인 내부에 다른 표고 제약이 하나도 없으면 그 안이 마루(바깥이 낮을 때) 또는 웅덩이(바깥이 높을 때)다. 바깥 사면의 국소 경사를 안쪽으로 연장하되 ±(간격−0.5m)로 제한한다 — 다음 등고선이 없다는 사실과 모순되지 않는 범위다. - 처리한 영역 수를 반환한다. + 처리한 영역의 마스크를 반환한다 — 이어지는 완화에서 함께 고정해야 한다. """ from scipy.ndimage import binary_dilation, binary_fill_holes, distance_transform_edt, label @@ -268,6 +286,7 @@ def _resolve_enclosed_interiors( band_cells = 8 limit = max(interval_m - 0.5, 0.5) resolved = 0 + handled = np.zeros(burned.shape, dtype=bool) for level in present: mask = burned == level interior = binary_fill_holes(mask) & ~mask @@ -301,10 +320,11 @@ def _resolve_enclosed_interiors( peak = float(inner.max()) offset = (limit / 2.0) * (inner[component] / peak) if peak > 0 else 0.0 surface[component] = level + direction * offset + handled |= component resolved += 1 if resolved: logger.info("도엽 서피스: 폐합 등고선 내부 %d곳을 사면 경사로 연장", resolved) - return resolved + return handled def build_sheet_surface_model( @@ -331,28 +351,25 @@ def build_sheet_surface_model( spec = grid_spec_from_bounds(x_min, y_min, x_max, y_max, SHEET_SURFACE_GRID_M) - from B04_PreProcess.B04_PreProcess_Engine_Watershed_Descent import rasterize_contours - # ① 2D — 등고 라인을 격자에 굽는다(라인 셀 = 표고, 그 외 NaN). - burned, levels = rasterize_contours(spec, features, None) - present = sorted({level for level in levels if bool((burned == level).any())}) + burned = _rasterize_contour_levels(spec, features) + present = sorted(np.unique(burned[np.isfinite(burned)]).tolist()) if len(present) < 2: logger.warning("도엽 서피스: 절취 범위 안에 등고선이 부족합니다.") return None - # 등고 간격(m) — 마루 캡 크기의 근거. 레벨이 하나뿐이면 5m(1:5,000 주곡선) 폴백. + # 등고 간격(m) — 마루 연장 상한의 근거. 레벨이 하나뿐이면 5m(주곡선) 폴백. interval_m = float(np.diff(np.array(present)).min()) if len(present) > 1 else 5.0 - # ② 2D — 계곡 구조선(하천중심선) 앵커 보간값을 같은 격자에 제약으로 굽는다. - stream_features = _load_features_metric(processed_dir, epsg, _STREAM_FILE) - if stream_features: - burned_stream = _burn_stream_anchors(spec, burned, stream_features) - logger.info("도엽 서피스: 계곡 구조선 제약 %d셀", burned_stream) - - # ③ 2D — 등고선 사이 거리 비례 보간(메시 없음). 여기서 나온 격자에서 1m 등고선을 + # ② 2D — 등고선 사이 거리 비례 보간(메시 없음). 여기서 나온 격자에서 1m 등고선을 # 뽑으므로 화면 등고선이 곧 2D 보간선이다(2026-08-30 사용자 지시). + # 제약은 **등고선만** 쓴다 — 계곡 구조선 앵커를 1m로 양자화해 섞었더니 제약 + # 표고가 29단→121단이 되어 계곡 주변만 1m 간격이 되고 보간선이 등간격을 + # 잃었다(2026-08-30 사용자 지적). V자 등고선이 계곡 하강을 이미 담고 있다. surface = _interpolate_between_contours(burned, spec.cell_m) - # ④ 폐합 등고선 안쪽(마루·웅덩이)을 바깥 사면 경사로 연장 (2026-08-30 사용자 확정). - _resolve_enclosed_interiors(burned, present, surface, spec.cell_m, interval_m) + # ③ 폐합 등고선 안쪽(마루·웅덩이)을 바깥 사면 경사로 연장 (2026-08-30 사용자 확정). + summits = _resolve_enclosed_interiors(burned, present, surface, spec.cell_m, interval_m) + # ④ 등고 라인(과 마루)을 고정한 채 완화 — 중심축 뭉침을 풀어 간격을 고르게 한다. + _relax_surface(surface, np.isfinite(burned) | summits, SHEET_SURFACE_RELAX_ITERATIONS) # DtmGridSampler 규약에 맞춰 y 오름차순으로 뒤집어 저장한다. x_coords = spec.cell_centers_x() diff --git a/config/config_system.py b/config/config_system.py index 56bb9ba2..64ed7dac 100644 --- a/config/config_system.py +++ b/config/config_system.py @@ -188,12 +188,10 @@ 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")) -# 가상 등고라인 반복 이분 횟수. 2회 = 5m 주곡선 → 2.5m → 1.25m 가상 등고 생성 -# (2026-08-30 사용자 지시 — 2D 사이 등고선을 먼저 만들고 3D화). -SHEET_SURFACE_MIDLINE_ROUNDS = int(os.getenv("SHEET_SURFACE_MIDLINE_ROUNDS", "2")) +# 등고선 고정 다듬기 반복 횟수. 거리 보간이 능선·계곡 중심축에 남기는 각진 자국만 +# 지우는 용도라 **적을수록 좋다** — 수렴시키면 harmonic 해가 되어 마루가 눌린다 +# (합성 원뿔 검증: 5회가 최적, 20회 이상 악화). 0이면 다듬지 않는다. +SHEET_SURFACE_RELAX_ITERATIONS = int(os.getenv("SHEET_SURFACE_RELAX_ITERATIONS", "5")) # 일반 사용자 WF1 자동 확정 기본값 SURFACE_CONFIRM_DEFAULT_FILTER = os.getenv("SURFACE_CONFIRM_DEFAULT_FILTER", "csf")