From c0f4c766b64b6dee1b0377bc0c4cfe20daf41bd4 Mon Sep 17 00:00:00 2001 From: umsangdon Date: Tue, 3 Nov 2026 13:46:33 +0900 Subject: [PATCH] =?UTF-8?q?feat(762-1):=20=EC=A7=80=EB=A9=B4=EC=A0=90=20?= =?UTF-8?q?=EB=8B=A4=EB=93=AC=EA=B8=B0(SMRF=20=EA=BC=B4)=20=C2=B7=20?= =?UTF-8?q?=EA=B3=A0=EB=A5=BC=20=EA=B0=88=EB=9E=98=20=EC=97=85=EC=B2=B4=20?= =?UTF-8?q?=EB=B6=84=EB=A5=98=20+=20=EB=8B=A4=EB=93=AC=EA=B8=B0=20=C2=B7?= =?UTF-8?q?=20=EC=A0=90=EC=A7=84=20TIN=20+=20=EB=8B=A4=EB=93=AC=EA=B8=B0?= =?UTF-8?q?=20=C2=B7=20=EC=A0=9C=EC=95=88=20=ED=9D=90=EB=A6=84=20=EC=8A=A4?= =?UTF-8?q?=EC=9C=84=EC=B9=98(=EA=B8=B0=EB=B3=B8=20=EB=81=94)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 새 Engine_Filter_Refine — 1 m 칸 최저 · 낮은 잡점(뒤집어 열림) · 창 키움 열림(최대 18 m · 기울기 0.3) · 걸린 칸 빼고 임시 지면 다시 그림 · 점마다 0.5 m + 1.25 × 기울기 - 필터 이름 classification_refine · ptd_refine — Ground _FILTERS · Schema · sampler SurfaceFilter 에 더함 · 엔진은 바탕 마스크 저장본을 다듬음(점진 TIN 두 번 안 돌림) - 스위치 SURFACE_REFINE_FLOW(프로세스 환경 변수 · 기본 끔) — 켜면 자동 전처리 · 확정 기본 = 다듬기 판 · 끄면 지금 길 그대로 - 기울기 0.3 = 용화 훑기에서 고름(PDAL 0.15 는 급경사 능선을 더 깎음 · 0.45 · 0.6 은 덤불 혹 더 남김) Co-Authored-By: Claude Opus 5.5 --- B04_PreProcess/B04_PreProcess_Engine.py | 32 +++-- .../B04_PreProcess_Engine_Filter_Refine.py | 116 ++++++++++++++++++ .../B04_PreProcess_Engine_Ground.py | 39 +++++- B04_PreProcess/B04_PreProcess_Schema.py | 7 +- common_util/common_util_surface_sampler.py | 2 + config/config_system_terrain.py | 14 ++- 6 files changed, 189 insertions(+), 21 deletions(-) create mode 100644 B04_PreProcess/B04_PreProcess_Engine_Filter_Refine.py diff --git a/B04_PreProcess/B04_PreProcess_Engine.py b/B04_PreProcess/B04_PreProcess_Engine.py index becc04c5e..37f9b0126 100644 --- a/B04_PreProcess/B04_PreProcess_Engine.py +++ b/B04_PreProcess/B04_PreProcess_Engine.py @@ -14,6 +14,7 @@ from typing import Any import numpy as np from B04_PreProcess.B04_PreProcess_Engine_Ground import ( + base_filter, build_ground_masks, resolve_auto_source_filters, run_ground_filter, @@ -282,24 +283,29 @@ def run_surface_analysis( source_filters = resolve_auto_source_filters(data, processed_dir) logger.info("B04 자동 전처리 기본 필터: %s", ", ".join(source_filters)) masks: dict[str, np.ndarray] = {} - for filter_key in source_filters: + + def _mask_for(filter_key: str) -> np.ndarray: mask_path = processed_dir / f"mask_{filter_key}.npy" - mask: np.ndarray | None = None if not rebuild and mask_path.is_file(): stored_mask = np.load(mask_path) if len(stored_mask) == total_points: - mask = np.asarray(stored_mask, dtype=bool) logger.info("B04 지면 필터 캐시 재사용: %s", mask_path.name) - if mask is None: - step_started = time.monotonic() - mask = np.asarray(run_ground_filter(filter_key, data), dtype=bool) - np.save(mask_path, mask) - logger.info( - "B04 지면 필터 계산 완료: %s (%.1fs)", - filter_key, - time.monotonic() - step_started, - ) - masks[filter_key] = mask + return np.asarray(stored_mask, dtype=bool) + # 다듬기 판(762)은 바탕 필터 마스크(저장본)를 다듬음 — 점진 TIN 을 두 번 돌리지 않음 + base = base_filter(filter_key) + base_mask = (masks.get(base) if base in masks else _mask_for(base)) if base else None + step_started = time.monotonic() + mask = np.asarray(run_ground_filter(filter_key, data, base_mask), dtype=bool) + np.save(mask_path, mask) + logger.info( + "B04 지면 필터 계산 완료: %s (%.1fs)", + filter_key, + time.monotonic() - step_started, + ) + return mask + + for filter_key in source_filters: + masks[filter_key] = _mask_for(filter_key) ground_summary = summarize_masks(data, masks) # 지면점 비율은 필터가 조용히 실패해도 유일하게 드러나는 신호다 — 반드시 남긴다. for filter_key, entry in ground_summary.items(): diff --git a/B04_PreProcess/B04_PreProcess_Engine_Filter_Refine.py b/B04_PreProcess/B04_PreProcess_Engine_Filter_Refine.py new file mode 100644 index 000000000..048da42d8 --- /dev/null +++ b/B04_PreProcess/B04_PreProcess_Engine_Filter_Refine.py @@ -0,0 +1,116 @@ +"""B04 지면점 다듬기 — 필터가 고른 지면점에서 덤불 · 풀 혹과 낮은 잡점을 한 번 더 거름(762-1). + +SMRF 꼴(Pingel 2013 · PDAL filters.smrf) — 입력 = 필터 지면 마스크 · 출력 = 줄어든 마스크. +① 1 m 칸 최저면 → 뒤집어 열림(낮은 잡점 칸) +② 창 반지름 1 칸부터 최대까지 키우며 열림 · 문턱 = 기울기 × 칸 × 반지름(물체 칸) +③ 걸린 칸 빼고 남은 칸으로 임시 지면 다시 그림(줄여 메우고 늘려 붓는 다층 보간) +④ 점마다 |z − 임시 지면| ≤ 문턱 + 배수 × 그 자리 기울기 이면 지면. +""" + +from typing import Any + +import numpy as np +from scipy import ndimage + +from config.config_system import ( + SURFACE_REFINE_CELL_M, + SURFACE_REFINE_LOW_SLOPE, + SURFACE_REFINE_SCALAR, + SURFACE_REFINE_SLOPE, + SURFACE_REFINE_THRESHOLD_M, + SURFACE_REFINE_WINDOW_M, +) + +_RELAX_STEPS = 30 + + +def _fill_nearest(grid: np.ndarray, missing: np.ndarray) -> np.ndarray: + if not missing.any(): + return grid + index = ndimage.distance_transform_edt(missing, return_distances=False, return_indices=True) + return grid[tuple(index)] + + +def _progressive(surface: np.ndarray, cell: float, slope: float, max_radius: int) -> np.ndarray: + """창을 한 칸씩 키운 열림이 한 단계에 문턱(기울기 × 칸 × 반지름)보다 깎은 칸.""" + flagged = np.zeros(surface.shape, dtype=bool) + last = surface + for radius in range(1, max_radius + 1): + # ponytail: 정사각 창(분리형 최소 · 최대 거름) — 원판보다 빠름 · 모서리 방향만 조금 넓음 + # (네 방향 선분 열림은 능선을 더 지키나 덤불 판을 3 배 남겨 버림 — 762 훑기) + current = ndimage.grey_opening(last, size=(2 * radius + 1, 2 * radius + 1)) + flagged |= (last - current) > slope * cell * radius + last = current + return flagged + + +def _inpaint(values: np.ndarray, known: np.ndarray) -> np.ndarray: + """빈 칸을 아는 칸으로 메움 — 2 배씩 줄여 메운 뒤 늘려 붓고 이웃 평균으로 폄.""" + if known.all(): + return values + if not known.any(): + return np.zeros_like(values) + rows, cols = values.shape + pad = ((0, rows % 2), (0, cols % 2)) + weight = np.pad(known.astype(np.float64), pad) + total = np.pad(np.where(known, values, 0.0), pad) + shape = (weight.shape[0] // 2, 2, weight.shape[1] // 2, 2) + coarse_weight = weight.reshape(shape).sum(axis=(1, 3)) + coarse_known = coarse_weight > 0 + coarse = total.reshape(shape).sum(axis=(1, 3)) / np.maximum(coarse_weight, 1.0) + coarse = _inpaint(coarse, coarse_known) + up = ndimage.zoom(coarse, 2, order=1, mode="nearest", grid_mode=True)[:rows, :cols] + result = np.where(known, values, up) + missing = ~known + for _ in range(_RELAX_STEPS): + padded = np.pad(result, 1, mode="edge") + mean = (padded[:-2, 1:-1] + padded[2:, 1:-1] + padded[1:-1, :-2] + padded[1:-1, 2:]) / 4 + result[missing] = mean[missing] + return result + + +def refine_ground_mask( + xyz: Any, + mask: np.ndarray, + *, + cell: float = SURFACE_REFINE_CELL_M, + slope: float = SURFACE_REFINE_SLOPE, + window_m: float = SURFACE_REFINE_WINDOW_M, + threshold: float = SURFACE_REFINE_THRESHOLD_M, + scalar: float = SURFACE_REFINE_SCALAR, + low_slope: float = SURFACE_REFINE_LOW_SLOPE, +) -> np.ndarray: + """필터 지면 마스크를 다듬어 돌려줌 — 지면이던 점만 지면으로 남음(늘지 않음).""" + mask = np.asarray(mask, dtype=bool) + index = np.flatnonzero(mask) + refined = np.zeros(mask.shape, dtype=bool) + if index.size < 3: + refined[index] = True + return refined + points = np.asarray(xyz[index], dtype=np.float64) + fx = (points[:, 0] - points[:, 0].min()) / cell + fy = (points[:, 1] - points[:, 1].min()) / cell + gx = fx.astype(np.int64) + gy = fy.astype(np.int64) + shape = (int(gy.max()) + 1, int(gx.max()) + 1) + zmin = np.full(shape, np.inf) + np.minimum.at(zmin, (gy, gx), points[:, 2]) + empty = ~np.isfinite(zmin) + filled = _fill_nearest(np.where(empty, 0.0, zmin), empty) + + # ① 낮은 잡점 — 뒤집은 면에서 한 칸 창 열림이 문턱 넘게 깎는 칸(PDAL 기본 기울기 5) + low = _progressive(-filled, cell, low_slope, 1) & ~empty + surface = _fill_nearest(filled, empty | low) + # ② 물체 — 창을 최대 반지름까지 키우며 + objects = _progressive(surface, cell, slope, max(1, int(round(window_m / cell)))) + # ③ 임시 지면 · 그 자리 기울기(rise / run) + provisional = _inpaint(filled, ~empty & ~low & ~objects) + grad_y, grad_x = np.gradient(provisional, cell) + gradient = np.hypot(grad_x, grad_y) + # ④ 점마다 판정 — 칸 가운데 기준 쌍선형 + coords = [fy - 0.5, fx - 0.5] + ground_z = ndimage.map_coordinates(provisional, coords, order=1, mode="nearest") + local_slope = ndimage.map_coordinates(gradient, coords, order=1, mode="nearest") + keep = np.abs(points[:, 2] - ground_z) <= threshold + scalar * local_slope + refined[index[keep]] = True + return refined diff --git a/B04_PreProcess/B04_PreProcess_Engine_Ground.py b/B04_PreProcess/B04_PreProcess_Engine_Ground.py index 790805f5a..76393f353 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_Ground.py +++ b/B04_PreProcess/B04_PreProcess_Engine_Ground.py @@ -1,6 +1,7 @@ """B04 지면 필터 오케스트레이션. -구조화된 포인트클라우드(structured.npz)에 대해 grid_min_z/csf/pmf/ransac +구조화된 포인트클라우드(structured.npz)에 대해 grid_min_z/csf/pmf/ransac/ptd/classification +(+ `_refine` 다듬기 판) 필터를 실행하여 지면 마스크 딕셔너리를 만든다. 필터 선택은 config의 source_filters를 따른다. """ @@ -20,7 +21,12 @@ from B04_PreProcess.B04_PreProcess_Engine_Filter_Grid import filter_grid_min_z from B04_PreProcess.B04_PreProcess_Engine_Filter_PMF import filter_pmf from B04_PreProcess.B04_PreProcess_Engine_Filter_PTD import filter_ptd from B04_PreProcess.B04_PreProcess_Engine_Filter_RANSAC import filter_ransac -from config.config_system import SURFACE_AUTO_FALLBACK_FILTER +from B04_PreProcess.B04_PreProcess_Engine_Filter_Refine import refine_ground_mask +from config.config_system import SURFACE_AUTO_FALLBACK_FILTER, SURFACE_REFINE_FLOW + +# 다듬기 판(762) — `{기본 필터}_refine` = 기본 필터 지면점을 SMRF 꼴로 한 번 더 거름 +REFINE_SUFFIX = "_refine" +_REFINABLE = ("classification", "ptd") # 필터 키 → 함수 매핑 _FILTERS = { @@ -30,10 +36,17 @@ _FILTERS = { "ransac": filter_ransac, "ptd": filter_ptd, "classification": filter_classification, + "classification_refine": lambda data: refine_ground_mask( + data["xyz"], filter_classification(data) + ), + "ptd_refine": lambda data: refine_ground_mask(data["xyz"], filter_ptd(data)), } # LAS에 지면분류가 실제로 들어 있을 때만 쓰는 필터 (미분류 LAS에서는 마스크가 빈다) -_CONDITIONAL_FILTERS = {"classification": has_classified_ground} +_CONDITIONAL_FILTERS = { + "classification": has_classified_ground, + "classification_refine": has_classified_ground, +} def available_filters() -> tuple[str, ...]: @@ -61,15 +74,29 @@ def resolve_auto_source_filters( usable = usable_conditional_filters(structured_data) if usable and not has_classified_coverage(structured_data, processed_dir): usable = [] - return [*usable[:1], SURFACE_AUTO_FALLBACK_FILTER] + filters = [*usable[:1], SURFACE_AUTO_FALLBACK_FILTER] + # 762 제안 흐름(스위치 켬) — 같은 기본 필터에 다듬기를 붙인 판을 지음 + if SURFACE_REFINE_FLOW: + filters = [f"{key}{REFINE_SUFFIX}" if key in _REFINABLE else key for key in filters] + return filters + + +def base_filter(filter_key: str) -> str | None: + """다듬기 판이면 바탕 필터 이름 — 엔진이 바탕 마스크 저장본을 다시 쓰게.""" + base = filter_key.removesuffix(REFINE_SUFFIX) + return base if base != filter_key and base in _REFINABLE else None def run_ground_filter( - filter_key: str, structured_data: dict[str, Any] | np.lib.npyio.NpzFile + filter_key: str, + structured_data: dict[str, Any] | np.lib.npyio.NpzFile, + base_mask: np.ndarray | None = None, ) -> np.ndarray: - """단일 지면 필터를 실행해 불리언 마스크를 반환한다.""" + """단일 지면 필터를 실행해 불리언 마스크를 반환한다 — 다듬기 판은 `base_mask` 를 다듬음.""" if filter_key not in _FILTERS: raise ValueError(f"알 수 없는 지면 필터입니다: {filter_key}") + if base_mask is not None and base_filter(filter_key): + return refine_ground_mask(structured_data["xyz"], base_mask) return np.asarray(_FILTERS[filter_key](structured_data), dtype=bool) diff --git a/B04_PreProcess/B04_PreProcess_Schema.py b/B04_PreProcess/B04_PreProcess_Schema.py index b13522af4..9225853d1 100644 --- a/B04_PreProcess/B04_PreProcess_Schema.py +++ b/B04_PreProcess/B04_PreProcess_Schema.py @@ -11,7 +11,12 @@ from config.config_system import ( # classification은 LAS가 지면분류를 싣고 있을 때 엔진이 자동으로 붙이지만, # 재계산 요청에서 이름으로 지정하는 것도 허용한다. -_ALLOWED_FILTERS = set(SURFACE_MODEL_SOURCE_FILTERS) | {"ransac", "classification"} +_ALLOWED_FILTERS = set(SURFACE_MODEL_SOURCE_FILTERS) | { + "ransac", + "classification", + "classification_refine", + "ptd_refine", +} _ALLOWED_METHODS = set(SURFACE_MODEL_PRECOMPUTE) diff --git a/common_util/common_util_surface_sampler.py b/common_util/common_util_surface_sampler.py index cc26565ac..eab838374 100644 --- a/common_util/common_util_surface_sampler.py +++ b/common_util/common_util_surface_sampler.py @@ -33,6 +33,8 @@ SurfaceFilter = Literal[ "pmf", "ptd", "ransac", + "classification_refine", + "ptd_refine", "sheet_tin_sheet", "sheet_tin", "sheet_biharmonic", diff --git a/config/config_system_terrain.py b/config/config_system_terrain.py index 104a7e0eb..8e0b260ee 100644 --- a/config/config_system_terrain.py +++ b/config/config_system_terrain.py @@ -94,6 +94,16 @@ SURFACE_MODEL_SOURCE_FILTERS = tuple( SURFACE_AUTO_METHODS = tuple(os.getenv("SURFACE_AUTO_METHODS", "dtm").split(",")) # 기본 필터는 입력 LAS를 보고 정한다 — 지면분류(class 2)가 있으면 그것을, 없으면 아래 값. SURFACE_AUTO_FALLBACK_FILTER = os.getenv("SURFACE_AUTO_FALLBACK_FILTER", "ptd") +# 762 제안 흐름 스위치 — 켜면 자동 전처리 · 확정 기본 = 기본 필터 + 다듬기(`_refine`) · 기본 끔. +SURFACE_REFINE_FLOW = os.getenv("SURFACE_REFINE_FLOW", "False").lower() == "true" +# 다듬기(SMRF 꼴) — 칸 · 물체 기울기(rise/run · PDAL 0.15 는 급경사 능선을 깎아 0.3) · 최대 창 · +# 판정 문턱 · 기울기 배수 · 낮은 잡점 기울기 — 나머지는 PDAL 기본. +SURFACE_REFINE_CELL_M = float(os.getenv("SURFACE_REFINE_CELL_M", "1.0")) +SURFACE_REFINE_SLOPE = float(os.getenv("SURFACE_REFINE_SLOPE", "0.3")) +SURFACE_REFINE_WINDOW_M = float(os.getenv("SURFACE_REFINE_WINDOW_M", "18.0")) +SURFACE_REFINE_THRESHOLD_M = float(os.getenv("SURFACE_REFINE_THRESHOLD_M", "0.5")) +SURFACE_REFINE_SCALAR = float(os.getenv("SURFACE_REFINE_SCALAR", "1.25")) +SURFACE_REFINE_LOW_SLOPE = float(os.getenv("SURFACE_REFINE_LOW_SLOPE", "5.0")) # dtm을 먼저 빌드해야 TIN 등고선 사전 캐시가 dtm footprint를 참조할 수 있다 (PLAN D-6) SURFACE_MODEL_PRECOMPUTE = tuple( os.getenv("SURFACE_MODEL_PRECOMPUTE", "dtm,tin,nurbs,implicit,meshfree").split(",") @@ -177,7 +187,9 @@ SHEET_SURFACE_AUTO_METHODS = [ ] # 일반 사용자 WF1 자동 확정 기본값 -SURFACE_CONFIRM_DEFAULT_FILTER = os.getenv("SURFACE_CONFIRM_DEFAULT_FILTER", "ptd") +SURFACE_CONFIRM_DEFAULT_FILTER = os.getenv( + "SURFACE_CONFIRM_DEFAULT_FILTER", "ptd_refine" if SURFACE_REFINE_FLOW else "ptd" +) # 확정 표현 = implicit · 평활 끔 — WF1 자동 짓기도 이 표현을 기본 필터에 지음(요구 558). SURFACE_CONFIRM_DEFAULT_METHOD = os.getenv("SURFACE_CONFIRM_DEFAULT_METHOD", "implicit") SURFACE_CONFIRM_DEFAULT_SMOOTH = (