From 8fedbd77acc40933871f552b36e4fe93dc8ce35f Mon Sep 17 00:00:00 2001 From: umsangdon Date: Sun, 30 Aug 2026 16:31:41 +0900 Subject: [PATCH] =?UTF-8?q?feat(B04):=20=EB=AF=B8=EC=B0=A9=EC=88=98?= =?UTF-8?q?=EC=98=80=EB=8D=98=20ANUDEM=ED=98=95=C2=B7=EB=8B=A4=EC=A4=91?= =?UTF-8?q?=ED=95=B4=EC=83=81=EB=8F=84=20=EB=B3=B4=EA=B0=84=EC=9D=84=20?= =?UTF-8?q?=EB=8D=94=ED=95=B4=208=EC=A2=85=EC=9C=BC=EB=A1=9C=20=EB=A7=8C?= =?UTF-8?q?=EB=93=A0=EB=8B=A4?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - ANUDEM형: 등고선의 국소 최대 곡률점(코너)에서 능선·계곡 구조선을 뽑아 (굽은 안쪽이 높으면 계곡, 낮으면 능선) 이웃 표고끼리 이어 앵커로 심고, 박판 스플라인으로 풀고, 형태학적 재구성으로 가짜 웅덩이를 메운다. Topo to Raster가 밟는 세 단계를 그대로 옮긴 것. 메운 웅덩이는 평탄해지지 않게 안쪽으로 cm 단위 경사를 남긴다(Garbrecht·Martz 평탄면 해소). - 다중해상도: 1/8→1/4→1/2→1 격자로 내려오며 앞 단계 해를 시작값으로 쓴다. - 박판 솔버를 solve_min_curvature()로 분리해 두 방식이 함께 쓴다. - 화면 버튼 순서에 두 방식을 넣었다(버튼이 곧 하이드·쇼). 실측(c1bb453f, 8종 83s) 평탄 셀: TPS 0.003 / 다중해상도 0.003 / 거리비례 0.014 / TIN 곡면 0.026 / ANUDEM형 0.046 / 라플라스 0.235 / TIN 선형 10.36 / IDW 69.13 %. LiDAR 대비 |dz|는 2.59~3.00m로 차이가 작다. Co-Authored-By: Claude Opus 5 (1M context) --- .../B04_PreProcess_Engine_SheetMethods.py | 221 +++++++++++++++++- B04_PreProcess/B04_PreProcess_UI_Page.ts | 2 + config/config_system.py | 3 +- 3 files changed, 218 insertions(+), 8 deletions(-) diff --git a/B04_PreProcess/B04_PreProcess_Engine_SheetMethods.py b/B04_PreProcess/B04_PreProcess_Engine_SheetMethods.py index b01bf7ac..10d11763 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_SheetMethods.py +++ b/B04_PreProcess/B04_PreProcess_Engine_SheetMethods.py @@ -10,6 +10,7 @@ """ import logging +import warnings from typing import Any, Callable import numpy as np @@ -20,6 +21,8 @@ logger = logging.getLogger(__name__) SHEET_METHOD_LABELS: dict[str, str] = { "distance": "거리비례", "biharmonic": "TPS(박판)", + "anudem": "ANUDEM형", + "multires": "다중해상도", "laplace": "라플라스", "tin": "TIN 선형", "clough": "TIN 곡면", @@ -129,16 +132,20 @@ def build_laplace(spec: Any, burned: np.ndarray, features: Any, cell_m: float) - # ── ③ 박판 스플라인(중조화) ────────────────────────────────────────────────── -def build_biharmonic(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray: - """등고선을 고정하고 Δ²z=0(박판 스플라인)을 켤레기울기법으로 푼다. +def solve_min_curvature(constrained: np.ndarray, guess: np.ndarray) -> np.ndarray: + """제약 셀을 고정하고 Δ²z=0(박판 스플라인)을 최소곡률 최소제곱으로 푼다. ANUDEM/Topo to Raster가 쓰는 박판 스플라인과 같은 연산자다. 라플라스와 달리 z=r을 그대로 통과시켜 능선·마루가 눌리지 않고, 경사가 등고선 너머로 자연스럽게 이어진다. + + `constrained`: 값이 고정된 셀(등고 라인, 필요하면 구조선 앵커) — 그 외는 NaN. + `guess`: 시작값이자 감쇠 기준(보통 거리 보간 결과). """ from scipy.sparse.linalg import LinearOperator, lsmr - guess = build_distance(spec, burned, features, cell_m).astype(np.float64) - # 등고 라인 + **격자 테두리**를 고정한다. 최외곽 등고선 바깥이 통째로 자유면 + guess = guess.astype(np.float64) + burned = constrained + # 제약 셀 + **격자 테두리**를 고정한다. 최외곽 등고선 바깥이 통째로 자유면 # 1차함수가 Δ²의 영공간에 남아 해가 하나로 정해지지 않고 켤레기울기가 발산한다 # (2026-08-30 실측: |Δz| 2092m). 테두리는 거리 보간값으로 묶는다. fixed = np.isfinite(burned) @@ -148,7 +155,8 @@ def build_biharmonic(spec: Any, burned: np.ndarray, features: Any, cell_m: float if not free.any(): return guess.astype(np.float32) index = np.flatnonzero(free.ravel()) - base = np.where(fixed, np.nan_to_num(guess), 0.0) + values = np.where(np.isfinite(burned), burned, guess) + base = np.where(fixed, np.nan_to_num(values), 0.0) # Δ²z=0을 정규방정식(CG)으로 풀면 조건수가 격자변 4제곱이라 발산한다(실측). # 대신 **최소곡률** 최소제곱으로 세운다 — 자유 셀에 대해 ‖Δz‖를 최소화하며, @@ -185,7 +193,7 @@ def build_biharmonic(spec: Any, burned: np.ndarray, features: Any, cell_m: float solution, info = result[0], result[1] surface = base.copy() surface.ravel()[index] = solution - surface[fixed] = guess[fixed] + surface[fixed] = values[fixed] # 안전장치 — 발산하면 조용히 틀린 지형을 넘기지 말고 거리 보간으로 되돌린다. drift = float(np.nanmax(np.abs(surface - guess))) span = float(np.nanmax(guess) - np.nanmin(guess)) @@ -196,10 +204,15 @@ def build_biharmonic(spec: Any, burned: np.ndarray, features: Any, cell_m: float info, ) return guess.astype(np.float32) - logger.info("도엽 서피스(TPS): 켤레기울기 수렴 info=%s, 최대 변화 %.2fm", info, drift) + logger.info("도엽 서피스(TPS): 최소제곱 info=%s, 최대 변화 %.2fm", info, drift) return surface.astype(np.float32) +def build_biharmonic(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray: + """등고선만 제약으로 둔 박판 스플라인.""" + return solve_min_curvature(burned, build_distance(spec, burned, features, cell_m)) + + # ── ④ TIN 선형 / ⑤ TIN 곡면 ───────────────────────────────────────────────── def _triangulated( spec: Any, burned: np.ndarray, cell_m: float, smooth: bool, max_points: int = 120_000 @@ -266,9 +279,203 @@ def build_idw(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np return surface.reshape(burned.shape).astype(np.float32) +# ── ⑦ ANUDEM형 (구조선 + 배수 강제) ───────────────────────────────────────── +def _contour_corner_anchors( + spec: Any, burned: np.ndarray, guess: np.ndarray, cell_m: float +) -> tuple[np.ndarray, np.ndarray] | None: + """등고선의 국소 최대 곡률점(코너)에서 능선·계곡 구조선 앵커를 만든다. + + ANUDEM은 등고선 자체의 곡률에서 능선·계곡망을 먼저 뽑아 흐름 구조를 세운다 + (Hutchinson 1988/89). 여기서도 같은 순서를 따른다. + + ① 라인마다 정점 곡률을 재 국소 최대점(V자 꼭짓점)을 고른다 + ② 굽은 안쪽이 더 높으면 **계곡**(등고선 V가 상류를 가리킴), 낮으면 **능선** + ③ 같은 종류의 코너를 이웃 표고끼리 이어 그 사이를 선형 보간해 앵커로 심는다 + + 앵커는 박판 해의 제약으로 들어가 계곡 바닥이 이어져 내려가고 능선 마루가 선다. + """ + from scipy.ndimage import label + from scipy.spatial import cKDTree + + levels = np.unique(burned[np.isfinite(burned)]) + if len(levels) < 2: + return None + interval = float(np.diff(levels).min()) + xs = spec.cell_centers_x() + ys = spec.cell_centers_y() + + corners: list[tuple[float, float, float, int]] = [] # x, y, level, +1 계곡 / -1 능선 + span = 6 # 곡률을 재는 정점 간격(px) — 짧으면 노이즈, 길면 꼭짓점을 놓친다 + for level in levels: + labelled, count = label(burned == level) + for component_id in range(1, count + 1): + line_rows, line_cols = np.nonzero(labelled == component_id) + if len(line_rows) < 3 * span: + continue + # 라인 셀을 한 줄로 세운다 — 좌표 정렬로 근사한다(정밀 추적은 과하다). + order = np.argsort(line_cols + line_rows * 1e-3) + path = np.column_stack([line_cols[order], line_rows[order]]).astype(np.float64) + before = np.roll(path, span, axis=0) + after = np.roll(path, -span, axis=0) + first = path - before + second = after - path + first_len = np.hypot(first[:, 0], first[:, 1]) + second_len = np.hypot(second[:, 0], second[:, 1]) + valid = (first_len > 1e-6) & (second_len > 1e-6) + cosine = np.ones(len(path)) + cosine[valid] = (first[valid] * second[valid]).sum(axis=1) / ( + first_len[valid] * second_len[valid] + ) + sharp = np.flatnonzero(valid & (cosine < 0.3)) # 70도 이상 꺾인 자리 + for i in sharp[:: max(1, span)]: + # 굽은 안쪽 방향 = 두 변 단위벡터 합의 반대 + inward = -(first[i] / first_len[i] + second[i] / second_len[i]) + norm = float(np.hypot(inward[0], inward[1])) + if norm < 1e-6: + continue + probe = path[i] + inward / norm * 6.0 + probe_col = int(round(probe[0])) + probe_row = int(round(probe[1])) + if not (0 <= probe_row < guess.shape[0] and 0 <= probe_col < guess.shape[1]): + continue + inside = float(guess[probe_row, probe_col]) + if not np.isfinite(inside) or abs(inside - level) < interval * 0.15: + continue + corners.append( + ( + float(xs[int(path[i, 0])]), + float(ys[int(path[i, 1])]), + float(level), + 1 if inside > level else -1, + ) + ) + if len(corners) < 4: + logger.info("도엽 서피스(ANUDEM형): 등고선 코너가 부족해 구조선을 건너뜁니다.") + return None + + array = np.asarray(corners, dtype=np.float64) + anchor_xy: list[np.ndarray] = [] + anchor_z: list[np.ndarray] = [] + reach = interval * 20.0 # 이보다 먼 코너는 같은 구조선으로 보지 않는다 + for level in levels[:-1]: + upper = level + interval + lower_set = array[np.abs(array[:, 2] - level) < 1e-6] + upper_set = array[np.abs(array[:, 2] - upper) < 1e-6] + if not len(lower_set) or not len(upper_set): + continue + tree = cKDTree(upper_set[:, :2]) + distance, index = tree.query(lower_set[:, :2], k=1) + for i in range(len(lower_set)): + j = int(index[i]) + if distance[i] > reach or lower_set[i, 3] != upper_set[j, 3]: + continue + start = lower_set[i, :2] + end = upper_set[j, :2] + steps = max(2, int(distance[i] / max(cell_m, 1e-6) / 4)) + fraction = np.linspace(0.0, 1.0, steps + 1)[1:-1] + if not len(fraction): + continue + anchor_xy.append(start + (end - start) * fraction[:, None]) + anchor_z.append(level + interval * fraction) + if not anchor_xy: + return None + return np.vstack(anchor_xy), np.concatenate(anchor_z) + + +def _enforce_drainage(surface: np.ndarray, epsilon: float = 0.01) -> int: + """가짜 웅덩이를 메운다 — ANUDEM의 배수 강제와 같은 목적. + + 등고선만으로 만든 면에는 흐름이 끊기는 웅덩이가 남는다. 형태학적 재구성(erosion) + 으로 채우되 완전 평탄해지지 않게 아주 작은 값을 얹는다. 채운 셀 수를 반환한다. + """ + from skimage.morphology import reconstruction + + if not np.isfinite(surface).all(): + return 0 + seed = np.full(surface.shape, float(surface.max()), dtype=np.float64) + seed[0, :] = surface[0, :] + seed[-1, :] = surface[-1, :] + seed[:, 0] = surface[:, 0] + seed[:, -1] = surface[:, -1] + filled = reconstruction(seed, surface.astype(np.float64), method="erosion") + raised = filled > surface + 1e-6 + if not raised.any(): + return 0 + # 그냥 채우면 웅덩이가 통째로 평탄해져 흐름 방향이 없어진다. 채운 영역 안쪽으로 + # 갈수록 아주 조금 높아지게 해서 물이 가장자리(넘침점)로 빠져나가게 둔다 + # (Garbrecht·Martz의 평탄면 해소를 간단히 옮긴 것 — 표고 변화는 cm 단위다). + from scipy.ndimage import distance_transform_edt + + inner = distance_transform_edt(raised) + surface[raised] = (filled[raised] + epsilon * inner[raised]).astype(surface.dtype) + return int(raised.sum()) + + +def build_anudem(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray: + """ANUDEM형 — 등고선 곡률에서 능선·계곡 구조선을 뽑아 제약에 더하고, 박판으로 풀고, + 가짜 웅덩이를 메운다. Topo to Raster가 밟는 세 단계를 그대로 옮긴 것이다.""" + guess = build_distance(spec, burned, features, cell_m).astype(np.float64) + constrained = burned.astype(np.float64).copy() + anchors = _contour_corner_anchors(spec, burned, guess, cell_m) + if anchors is not None: + xy, z = anchors + 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(constrained[row, col]) + constrained[row[free], col[free]] = z[free] + logger.info("도엽 서피스(ANUDEM형): 구조선 앵커 %d셀", int(free.sum())) + surface = solve_min_curvature(constrained, guess).astype(np.float32) + logger.info("도엽 서피스(ANUDEM형): 가짜 웅덩이 %d셀 메움", _enforce_drainage(surface)) + return surface + + +# ── ⑧ 다중해상도 (coarse → fine) ──────────────────────────────────────────── +def build_multires(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray: + """성긴 격자에서 풀고 점차 세밀화한다 — ANUDEM의 다중해상도 전략. + + 전체 형상은 성긴 격자에서 싸게 잡고, 세밀한 격자에서는 등고선 근처만 다듬는다. + 한 해상도에서만 풀 때보다 넓은 밴드가 고르게 퍼지고 값싸게 수렴한다. + """ + from scipy.ndimage import zoom + + surface: np.ndarray | None = None + for factor in (8, 4, 2, 1): + if factor == 1: + coarse = burned + else: + # 성긴 격자의 제약 — 블록 안 등고 라인의 평균 표고를 대표로 쓴다. + rows = burned.shape[0] // factor * factor + cols = burned.shape[1] // factor * factor + blocks = burned[:rows, :cols].reshape(rows // factor, factor, cols // factor, factor) + # 라인이 하나도 없는 블록은 NaN이 정상이라 경고를 삼킨다. + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + coarse = np.nanmean(blocks, axis=(1, 3)).astype(np.float32) + level = build_distance(spec, coarse, features, cell_m * factor) + if surface is not None: + # 앞 단계 해를 지금 해상도로 올려 절반씩 섞는다 — 성긴 단계의 넓은 추세를 + # 이어받되 이번 해상도의 등고선 정보를 덮지 않는다. + scale = (level.shape[0] / surface.shape[0], level.shape[1] / surface.shape[1]) + upscaled = zoom(surface, scale, order=1) + free = ~np.isfinite(coarse) + level[free] = (level[free] + upscaled[free]) * 0.5 + relax_laplace(level, np.isfinite(coarse), 8) + surface = level + assert surface is not None + if surface.shape != burned.shape: # 블록 자르기로 남은 가장자리 보정 + scale = (burned.shape[0] / surface.shape[0], burned.shape[1] / surface.shape[1]) + surface = zoom(surface, scale, order=1) + line = np.isfinite(burned) + surface[line] = burned[line] + return surface.astype(np.float32) + + SHEET_METHOD_BUILDERS: dict[str, Callable[[Any, np.ndarray, Any, float], np.ndarray]] = { "distance": build_distance, "biharmonic": build_biharmonic, + "anudem": build_anudem, + "multires": build_multires, "laplace": build_laplace, "tin": build_tin, "clough": build_clough, diff --git a/B04_PreProcess/B04_PreProcess_UI_Page.ts b/B04_PreProcess/B04_PreProcess_UI_Page.ts index 45ced88b..f6e366ac 100644 --- a/B04_PreProcess/B04_PreProcess_UI_Page.ts +++ b/B04_PreProcess/B04_PreProcess_UI_Page.ts @@ -47,6 +47,8 @@ const ROUTE_STAGE = ROUTES.B05_PROFILE; const SHEET_METHOD_ORDER = [ "distance", "biharmonic", + "anudem", + "multires", "laplace", "tin", "clough", diff --git a/config/config_system.py b/config/config_system.py index 6bda0157..06f61284 100644 --- a/config/config_system.py +++ b/config/config_system.py @@ -193,7 +193,8 @@ SHEET_SURFACE_GRID_M = float(os.getenv("SHEET_SURFACE_GRID_M", "1.0")) SHEET_SURFACE_METHODS = [ method.strip() for method in os.getenv( - "SHEET_SURFACE_METHODS", "distance,biharmonic,laplace,tin,clough,idw" + "SHEET_SURFACE_METHODS", + "distance,biharmonic,anudem,multires,laplace,tin,clough,idw", ).split(",") if method.strip() ]