refactor(767-1): 안 쓰는 지면 필터 · 서피스 지움 — 필터 넷(CSF · 격자 최저 · PMF · RANSAC) · 표현 셋(TIN + Taubin 평활 · NURBS · meshfree) · 수치지형도 비교 방식 다섯 · 안 부르는 길 정리
- 필터 = 업체 분류 · 점진 TIN · 둘의 다듬기 판만 · 분석 허용 목록 = Ground 등록부(available_filters) 한 벌 - 표현 = DTM · implicit 만 · 평활 = DTM 만 · config 기본 목록 · PRECOMPUTE · 설정 사전 맞춤(설정 지문 바뀜 — 768 에서 다시 지음) - 등고선 = 격자 · 굳힌 격자만(삼각망 굽기 · 외곽 맞춤 · implicit RBF 폴백 걷음) · 표본기 · 경로 찾기 tin · nurbs · meshfree 갈래 걷음 - 수치지형도 = 라플라스만(TIN 둘 · TPS · ANUDEM형 · 다중해상도 걷음) · B04 단추 하나 - 안 부르는 길 — B04 JSON 점군 GET · ground-stats GET + TS · B07 옛 계획평면도(라이다) 짓기 · lidar_source · LIDAR_ID · plan_lidar 종류 이름 · SOURCE_FILTER_LABELS · MODEL_METHODS(py) · MESH_* 설정 - 옛 manifest 의 지운 갈래 줄은 DB 로 되살리지 않음 · 옛 csf 행은 미리보기 · 등고선 그대로 열림 · 지운 이름 요청 = 400 - 낡은 글 — Ground · Pipeline 머리 · WF1 「없으면 csf」 · B07 음영 기본 「csf」 → 빈 값 · 등고선 폴백 「csf」 → 404 · B04 3D 보기 처음값 · B05 Schema · B06 Prebuild · 씨닝 설명 - 시험 지움 2 — test_b04_contour_raster · test_b04_contour_equivalence(지운 삼각망 굽기 · 없는 프로젝트를 봐 건너뛰던 것) - 시험 덜어냄 — test_b04_ground_filters CSF · 격자 최저 네 경우(지운 필터) - 시험 고침 — 허용 목록 · 형 검사에 걸린 csf · tin · grid_min_z · meshfree · tin_sheet 견본을 ptd · implicit · laplace 로(54_22 · 238_1 · 443 · 342_2 · 346_3 · 351_1 · 442_2 · 271 + helper · 344_1 · 133_1 · 343_2 · 637 · 754_1 · 763_1 깨진 tin → 깨진 implicit · 763_3 · b05 둘 · b06 · sheet_surface TIN(도엽선) 경우 지움) · 이름만 바꿈(273_2 · 343_1 · 343_3 · 624_1) Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5.5
@@ -1,57 +1,23 @@
|
||||
"""도엽등고선 → 표고 격자 보간 방식 모음 (비교용).
|
||||
"""도엽등고선 → 표고 격자 보간 — 라플라스(거리 비례 초기값 + 이웃 평균 다듬기).
|
||||
|
||||
문헌(Hutchinson 1988/89 ANUDEM; Chaplot 2006; Arun 2013 등)은 지형 복잡도·자료 밀도에
|
||||
따라 우열이 갈리며 단일 최적해가 없다고 본다. 그래서 방식을 하나로 고르지 않고 여기
|
||||
모아 두고 B04 화면에서 바꿔 가며 보게 한다(2026-08-30 사용자 지시).
|
||||
비교용 다섯 방식(TIN 둘 · TPS · ANUDEM형 · 다중해상도)은 걷음(608) — 확정 · 자동 = 라플라스.
|
||||
|
||||
각 builder는 `(spec, burned, features, cell_m) -> (R, C) float32` 격자를 돌려준다.
|
||||
builder는 `(spec, burned, features, cell_m) -> (R, C) float32` 격자를 돌려준다.
|
||||
`burned`는 등고 라인이 구워진 격자(라인 셀 = 표고, 그 외 NaN)다. 폐합 링 안쪽 처리와
|
||||
프리뷰·저장은 호출측(`_SheetSurface`)이 방식과 무관하게 똑같이 해 준다.
|
||||
"""
|
||||
|
||||
import logging
|
||||
import warnings
|
||||
from typing import Any, Callable
|
||||
|
||||
import numpy as np
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
# TIN(도엽선)이 물어 오는 격자 밖 여유 — 테두리가 삼각망 밖으로 나가지 않을 만큼만.
|
||||
SHEET_TIN_CLIP_MARGIN_M = 100.0
|
||||
|
||||
# 화면 버튼에 쓰는 이름 — 키는 surface_models.generation_params.source_filter 접미사다.
|
||||
SHEET_METHOD_LABELS: dict[str, str] = {
|
||||
"tin_sheet": "TIN(도엽선)",
|
||||
"tin": "TIN 격자",
|
||||
"biharmonic": "TPS(박판)",
|
||||
"anudem": "ANUDEM형",
|
||||
"multires": "다중해상도",
|
||||
"laplace": "라플라스",
|
||||
}
|
||||
|
||||
|
||||
def _laplacian(values: np.ndarray) -> np.ndarray:
|
||||
padded = np.pad(values, 1, mode="edge")
|
||||
return (
|
||||
padded[:-2, 1:-1] + padded[2:, 1:-1] + padded[1:-1, :-2] + padded[1:-1, 2:] - 4.0 * values
|
||||
)
|
||||
|
||||
|
||||
def _contour_vertices(burned: np.ndarray, spec: Any) -> tuple[np.ndarray, np.ndarray]:
|
||||
"""등고 라인 셀을 (N,2) 세계좌표와 표고로 바꾼다."""
|
||||
rows, cols = np.nonzero(np.isfinite(burned))
|
||||
xs = spec.cell_centers_x()[cols]
|
||||
ys = spec.cell_centers_y()[rows]
|
||||
return np.column_stack([xs, ys]), burned[rows, cols].astype(np.float64)
|
||||
|
||||
|
||||
def _grid_points(spec: Any) -> tuple[np.ndarray, np.ndarray]:
|
||||
grid_x, grid_y = np.meshgrid(spec.cell_centers_x(), spec.cell_centers_y())
|
||||
return grid_x, grid_y
|
||||
|
||||
|
||||
# ── ① 거리 비례 (버튼에서는 뺐지만 다른 방식의 초기추정으로 계속 쓴다) ───────
|
||||
# ── 거리 비례 (라플라스 초기값) ─────────────────────────────────────────────
|
||||
def build_distance(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
|
||||
"""가장 가까운 서로 다른 표고 두 라인 사이를 거리 비례로 나눈다.
|
||||
|
||||
@@ -120,394 +86,17 @@ def relax_laplace(surface: np.ndarray, fixed: np.ndarray, iterations: int) -> No
|
||||
surface[colour] = neighbours[colour]
|
||||
|
||||
|
||||
# ── ② 라플라스(조화) ─────────────────────────────────────────────────────────
|
||||
# ── 라플라스(조화) ─────────────────────────────────────────────────────────────
|
||||
def build_laplace(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
|
||||
"""등고선을 경계값으로 두고 Δz=0을 푼다.
|
||||
|
||||
면이 매끈해지지만 z=r(원뿔·능선)은 조화함수가 아니라 마루가 눌린다. 비교 기준으로
|
||||
남겨 둔다 — ANUDEM이 라플라스 대신 박판 스플라인을 쓰는 이유를 눈으로 보기 위함.
|
||||
면이 매끈해지지만 z=r(원뿔·능선)은 조화함수가 아니라 마루가 조금 눌린다.
|
||||
"""
|
||||
surface = build_distance(spec, burned, features, cell_m)
|
||||
relax_laplace(surface, np.isfinite(burned), 400)
|
||||
return surface
|
||||
|
||||
|
||||
# ── ③ 박판 스플라인(중조화) ──────────────────────────────────────────────────
|
||||
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 = guess.astype(np.float64)
|
||||
burned = constrained
|
||||
# 제약 셀 + **격자 테두리**를 고정한다. 최외곽 등고선 바깥이 통째로 자유면
|
||||
# 1차함수가 Δ²의 영공간에 남아 해가 하나로 정해지지 않고 켤레기울기가 발산한다
|
||||
# (2026-08-30 실측: |Δz| 2092m). 테두리는 거리 보간값으로 묶는다.
|
||||
fixed = np.isfinite(burned)
|
||||
fixed[0, :] = fixed[-1, :] = True
|
||||
fixed[:, 0] = fixed[:, -1] = True
|
||||
free = ~fixed & np.isfinite(guess)
|
||||
if not free.any():
|
||||
return guess.astype(np.float32)
|
||||
index = np.flatnonzero(free.ravel())
|
||||
values = np.where(np.isfinite(burned), burned, guess)
|
||||
base = np.where(fixed, np.nan_to_num(values), 0.0)
|
||||
|
||||
# Δ²z=0을 정규방정식(CG)으로 풀면 조건수가 격자변 4제곱이라 발산한다(실측).
|
||||
# 대신 **최소곡률** 최소제곱으로 세운다 — 자유 셀에 대해 ‖Δz‖를 최소화하며,
|
||||
# 그 정상해가 곧 Δ²z=0이다. 조건수가 제곱으로 줄어 LSMR이 안정적으로 푼다
|
||||
# (Briggs 1974의 최소곡률 격자화와 같은 목적함수).
|
||||
# 최소곡률만으로는 제약(등고선)에서 먼 영역이 정해지지 않아 해가 폭주한다.
|
||||
# ANUDEM의 거칠기 벌점과 같은 취지로 감쇠항을 붙여 거리 보간값에 묶어 둔다:
|
||||
# minimize ‖Δz‖² + λ‖z − 거리보간‖²
|
||||
# λ가 작을수록 더 매끈하고 클수록 거리 보간에 가깝다.
|
||||
total_cells = base.size
|
||||
free_count = len(index)
|
||||
damping = np.float64(np.sqrt(0.02))
|
||||
anchor = guess.ravel()[index]
|
||||
|
||||
def forward(vector: np.ndarray) -> np.ndarray:
|
||||
# 반드시 **선형**이어야 한다 — 고정 셀 기여(base)를 여기서 더하면 아핀이 되어
|
||||
# LSMR의 전제가 깨지고 해가 폭주한다. base 몫은 우변으로만 넘긴다.
|
||||
scattered = np.zeros_like(base)
|
||||
scattered.ravel()[index] = vector
|
||||
return np.concatenate([_laplacian(scattered).ravel(), damping * vector])
|
||||
|
||||
def adjoint(vector: np.ndarray) -> np.ndarray:
|
||||
curvature = _laplacian(vector[:total_cells].reshape(base.shape)).ravel()[index]
|
||||
return curvature + damping * vector[total_cells:]
|
||||
|
||||
rhs = np.concatenate([-_laplacian(base).ravel(), damping * anchor])
|
||||
linear = LinearOperator(
|
||||
(total_cells + free_count, free_count),
|
||||
matvec=forward,
|
||||
rmatvec=adjoint,
|
||||
dtype=np.float64,
|
||||
)
|
||||
result = lsmr(linear, rhs, x0=anchor, maxiter=400, atol=1e-8, btol=1e-8)
|
||||
solution, info = result[0], result[1]
|
||||
surface = base.copy()
|
||||
surface.ravel()[index] = solution
|
||||
surface[fixed] = values[fixed]
|
||||
# 안전장치 — 발산하면 조용히 틀린 지형을 넘기지 말고 거리 보간으로 되돌린다.
|
||||
drift = float(np.nanmax(np.abs(surface - guess)))
|
||||
span = float(np.nanmax(guess) - np.nanmin(guess))
|
||||
if not np.isfinite(drift) or drift > max(span, 1.0):
|
||||
logger.warning(
|
||||
"도엽 서피스(TPS): 해가 발산해(최대 %.1fm) 거리 보간으로 되돌립니다 (info=%s).",
|
||||
drift,
|
||||
info,
|
||||
)
|
||||
return guess.astype(np.float32)
|
||||
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 ───────────────────────────────────────────────────────────────────
|
||||
def _triangulate(
|
||||
spec: Any, shape_: tuple[int, int], points: np.ndarray, values: np.ndarray
|
||||
) -> np.ndarray:
|
||||
"""정점 구름을 Delaunay 삼각망 선형 보간해 격자로 편다. 삼각망 밖은 NaN."""
|
||||
from scipy.interpolate import LinearNDInterpolator
|
||||
|
||||
if len(points) < 3:
|
||||
return np.full(shape_, np.nan, dtype=np.float32)
|
||||
interpolator = LinearNDInterpolator(points, values)
|
||||
grid_x, grid_y = _grid_points(spec)
|
||||
surface = np.empty(shape_, dtype=np.float32)
|
||||
chunk = max(1, int(4_000_000 // max(spec.n_cols, 1)))
|
||||
for start in range(0, spec.n_rows, chunk):
|
||||
stop = min(start + chunk, spec.n_rows)
|
||||
surface[start:stop] = interpolator(grid_x[start:stop], grid_y[start:stop]).astype(
|
||||
np.float32
|
||||
)
|
||||
return surface
|
||||
|
||||
|
||||
def build_tin(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
|
||||
"""격자에 구운 등고 라인 셀을 Delaunay 삼각망 선형 보간한다.
|
||||
|
||||
정점이 셀 중심에 맞춰져 있어 1m 계단이 삼각망에 그대로 실린다. 같은 표고 정점
|
||||
3개로 이루어진 평탄 삼각형이 굴곡부·마루에 계단을 만든다. 비교 기준으로 남긴다.
|
||||
"""
|
||||
points, values = _contour_vertices(burned, spec)
|
||||
if len(points) > 120_000: # 삼각망 비용은 정점 수에 비례한다
|
||||
step = int(np.ceil(len(points) / 120_000))
|
||||
points, values = points[::step], values[::step]
|
||||
return _triangulate(spec, burned.shape, points, values)
|
||||
|
||||
|
||||
def build_tin_sheet(spec: Any, burned: np.ndarray, features: Any, cell_m: float) -> np.ndarray:
|
||||
"""**원본 5m 도엽 등고선 정점**을 그대로 이은 고전 TIN (2026-08-30 사용자 지시).
|
||||
|
||||
`tin`은 격자에 구운 라인 셀(=1m 계단으로 뭉개진 정점)을 쓰지만, 이쪽은 벡터
|
||||
등고선의 정점을 좌표 그대로 쓴다 — 도면 등고선을 삼각망으로 잇는 측량 관행 그대로다.
|
||||
격자 밖 등고선은 물지 않는다(삼각망 비용만 커지고 결과는 같다). 다만 격자 테두리가
|
||||
삼각망 밖으로 나가지 않도록 여유를 두고 자른다.
|
||||
|
||||
길이 필터는 두지 않는다 — 짧은 봉우리 폐합 등고선을 버리면 마루가 통째로 평평해진다.
|
||||
"""
|
||||
from shapely.geometry import shape as to_shape
|
||||
|
||||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import (
|
||||
ELEVATION_KEYS,
|
||||
iter_linestrings,
|
||||
)
|
||||
|
||||
xs, ys = spec.cell_centers_x(), spec.cell_centers_y()
|
||||
margin = max(SHEET_TIN_CLIP_MARGIN_M, cell_m * 2.0)
|
||||
x_lo, x_hi = xs[0] - margin, xs[-1] + margin
|
||||
y_lo, y_hi = ys[-1] - margin, ys[0] + margin
|
||||
|
||||
coords: list[np.ndarray] = []
|
||||
levels: list[np.ndarray] = []
|
||||
for feature in features or []:
|
||||
properties = feature.get("properties") or {}
|
||||
elevation = next(
|
||||
(float(properties[key]) for key in ELEVATION_KEYS if properties.get(key) is not None),
|
||||
None,
|
||||
)
|
||||
geometry = feature.get("geometry")
|
||||
if elevation is None or not geometry:
|
||||
continue
|
||||
try:
|
||||
parsed = to_shape(geometry)
|
||||
except Exception: # noqa: BLE001 — 손상된 피처는 건너뛴다
|
||||
continue
|
||||
for line in iter_linestrings(parsed):
|
||||
point = np.asarray(line.coords, dtype=np.float64)[:, :2]
|
||||
inside = (
|
||||
(point[:, 0] >= x_lo)
|
||||
& (point[:, 0] <= x_hi)
|
||||
& (point[:, 1] >= y_lo)
|
||||
& (point[:, 1] <= y_hi)
|
||||
)
|
||||
if not inside.any():
|
||||
continue
|
||||
coords.append(point[inside])
|
||||
levels.append(np.full(int(inside.sum()), elevation, dtype=np.float64))
|
||||
|
||||
if not coords:
|
||||
logger.warning("도엽 서피스: TIN(도엽선)에 쓸 등고선 정점이 없습니다.")
|
||||
return np.full(burned.shape, np.nan, dtype=np.float32)
|
||||
|
||||
points = np.vstack(coords)
|
||||
values = np.concatenate(levels)
|
||||
# 도엽 이음매에서 같은 정점이 겹쳐 들어온다 — Qhull 비용만 늘어 미리 접는다.
|
||||
_, unique = np.unique(np.round(points, 3), axis=0, return_index=True)
|
||||
points, values = points[unique], values[unique]
|
||||
logger.info("도엽 서피스: TIN(도엽선) 정점 %d개", len(points))
|
||||
return _triangulate(spec, burned.shape, points, values)
|
||||
|
||||
|
||||
# ── ⑦ 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]] = {
|
||||
"tin_sheet": build_tin_sheet,
|
||||
"tin": build_tin,
|
||||
"biharmonic": build_biharmonic,
|
||||
"anudem": build_anudem,
|
||||
"multires": build_multires,
|
||||
"laplace": build_laplace,
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user