diff --git a/B04_PreProcess/B04_PreProcess_Engine.py b/B04_PreProcess/B04_PreProcess_Engine.py index 3789d2df..3f470d98 100644 --- a/B04_PreProcess/B04_PreProcess_Engine.py +++ b/B04_PreProcess/B04_PreProcess_Engine.py @@ -15,6 +15,7 @@ import numpy as np from B04_PreProcess.B04_PreProcess_Engine_Ground import ( build_ground_masks, + detect_extra_filters, run_ground_filter, summarize_masks, ) @@ -22,7 +23,7 @@ from B04_PreProcess.B04_PreProcess_Engine_Pipeline import build_all_terrain_mode from B04_PreProcess.B04_PreProcess_Engine_Structurize import structurize_las from common_util.common_util_atomic import atomic_write_npz from common_util.common_util_json import atomic_write_json -from config.config_system import build_surface_model_config +from config.config_system import SURFACE_GROUND_RATIO_WARN, build_surface_model_config logger = logging.getLogger(__name__) @@ -169,10 +170,15 @@ def run_surface_analysis( "y_min": float(bounds[1, 0]), "y_max": float(bounds[1, 1]), } - data = {"xyz": xyz, "bounds": bounds} + data = {"xyz": xyz, "bounds": bounds, "classification": structured["classification"]} # 2. 지면 필터 실행 — mask_{filter}.npy 영구 캐시 우선 재사용 (PLAN A-1) _report(40, "ground_filter", "지면 필터 적용 중") + # LAS가 지면분류를 싣고 왔으면 필터 하나로 더 쓴다 — 미분류 LAS면 조용히 빠진다. + extra_filters = detect_extra_filters(data, source_filters) + if extra_filters: + source_filters = list(source_filters) + extra_filters + logger.info("B04 LAS 자체 지면분류 감지 — 필터 추가: %s", ", ".join(extra_filters)) masks: dict[str, np.ndarray] = {} for filter_key in source_filters: mask_path = processed_dir / f"mask_{filter_key}.npy" @@ -193,6 +199,17 @@ def run_surface_analysis( ) masks[filter_key] = mask ground_summary = summarize_masks(data, masks) + # 지면점 비율은 필터가 조용히 실패해도 유일하게 드러나는 신호다 — 반드시 남긴다. + for filter_key, entry in ground_summary.items(): + ratio = float(entry["ground_ratio"]) + log = logger.warning if ratio < SURFACE_GROUND_RATIO_WARN else logger.info + log( + "B04 지면점 %s: %d / %d (%.2f%%)", + filter_key, + entry["ground_point_count"], + entry["total_point_count"], + ratio * 100, + ) for filter_key, mask in masks.items(): cache_path = processed_dir / f"ground_points_{filter_key}.npz" if not rebuild and cache_path.is_file(): diff --git a/B04_PreProcess/B04_PreProcess_Engine_Filter_CSF.py b/B04_PreProcess/B04_PreProcess_Engine_Filter_CSF.py index 3b861d77..8b77d808 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_Filter_CSF.py +++ b/B04_PreProcess/B04_PreProcess_Engine_Filter_CSF.py @@ -14,6 +14,7 @@ from config.config_system import ( SURFACE_CSF_CLASS_THRESHOLD_M, SURFACE_CSF_CLOTH_RESOLUTION_M, SURFACE_CSF_ITERATIONS, + SURFACE_CSF_MAX_ITERATIONS, SURFACE_CSF_RIGIDNESS, SURFACE_CSF_SLOPE_SMOOTH, SURFACE_CSF_SLOPE_SMOOTH_THRESHOLD_M, @@ -35,6 +36,7 @@ def filter_csf( iterations: int = SURFACE_CSF_ITERATIONS, slope_smooth: bool = SURFACE_CSF_SLOPE_SMOOTH, slope_smooth_threshold: float = SURFACE_CSF_SLOPE_SMOOTH_THRESHOLD_M, + max_iterations: int = SURFACE_CSF_MAX_ITERATIONS, ) -> np.ndarray: """CSF로 지면 포인트의 불리언 마스크를 반환한다.""" if not math.isfinite(cloth_resolution) or cloth_resolution <= 0: @@ -47,6 +49,8 @@ def filter_csf( raise ValueError("분류 임계값은 0 이상의 유한한 값이어야 합니다.") if iterations <= 0: raise ValueError("반복 횟수는 1 이상이어야 합니다.") + if max_iterations < iterations: + raise ValueError("반복 상한은 정착 여유 반복 수보다 작을 수 없습니다.") xyz = np.asarray(structured_data["xyz"], dtype=np.float64) if xyz.ndim != 2 or xyz.shape[1] != 3: @@ -78,9 +82,15 @@ def filter_csf( collision_grid[collision_grid == -np.inf] = 0.0 # 4. 천 시뮬레이션 반복 루프 (물리 하강) + # 반복 수는 현장 기복에서 뽑는다 — 고정 횟수는 곧 고정 하강 예산이라 + # (0.3185m × 150회 = 47.8m) 기복이 그보다 큰 산악지에서 천이 지면에 닿지 + # 못하고 공중에 멈춘다. 낙하 구간 + 스프링 정착 여유(iterations)로 잡는다. gravity = _GRAVITY_BASE * time_step spring_coeff = _RIGIDNESS_SPRING_COEFF[rigidness] - for _ in range(iterations): + required_drop = start_height - float(np.min(collision_grid)) + total_iterations = int(math.ceil(max(required_drop, 0.0) / gravity)) + iterations + total_iterations = min(total_iterations, max_iterations) + for _ in range(total_iterations): cloth_z -= gravity cloth_z = np.maximum(cloth_z, collision_grid) @@ -102,9 +112,11 @@ def filter_csf( height_diff = np.abs(inverted_zs - simulated_inverted_z) mask = height_diff <= class_threshold - # 6. 수목 노이즈 2차 필터 보정 + # 6. 수목 노이즈 2차 필터 보정 — 셀 지면 후보(반전고 최대 = 원표고 최저)보다 + # 얼마나 위에 떠 있는지를 본다. 반전 공간에서 높이는 (지면 후보 - 자기 반전고)다. + # 반대로 빼면 항상 0 이하가 나와 조건이 늘 참이 된다(= 필터가 무력화). if slope_smooth: - local_min_z = collision_grid[gy, gx] - mask = mask & ((inverted_zs - local_min_z) < slope_smooth_threshold) + local_ground_z = collision_grid[gy, gx] + mask = mask & ((local_ground_z - inverted_zs) < slope_smooth_threshold) return np.asarray(mask, dtype=bool) diff --git a/B04_PreProcess/B04_PreProcess_Engine_Filter_Classification.py b/B04_PreProcess/B04_PreProcess_Engine_Filter_Classification.py new file mode 100644 index 00000000..b3475e88 --- /dev/null +++ b/B04_PreProcess/B04_PreProcess_Engine_Filter_Classification.py @@ -0,0 +1,40 @@ +"""B04 LAS 자체 지면분류(ASPRS class 2) 필터. + +계산이 없다 — LAS 제작 업체가 이미 분류해 둔 지면점을 그대로 읽는다. +분류가 없는 LAS(classification 전부 0)에서는 마스크가 비므로, +쓸 수 있는지는 has_classified_ground()로 먼저 판정한다. +""" + +from typing import Any + +import numpy as np + +from config.config_system import SURFACE_CLASSIFIED_GROUND_MIN_RATIO + +# ASPRS 표준 지면 클래스 코드 +GROUND_CLASS_CODE = 2 + + +def filter_classification( + structured_data: dict[str, Any] | np.lib.npyio.NpzFile, +) -> np.ndarray: + """LAS classification이 지면(2)인 포인트의 불리언 마스크를 반환한다.""" + if "classification" not in structured_data: + raise ValueError("구조화 데이터에 classification 배열이 없습니다.") + classification = np.asarray(structured_data["classification"]) + return np.asarray(classification == GROUND_CLASS_CODE, dtype=bool) + + +def has_classified_ground( + structured_data: dict[str, Any] | np.lib.npyio.NpzFile, + minimum_ratio: float = SURFACE_CLASSIFIED_GROUND_MIN_RATIO, +) -> bool: + """LAS에 쓸 만한 지면분류가 들어 있는지 판정한다.""" + if "classification" not in structured_data: + return False + classification = np.asarray(structured_data["classification"]) + total = int(classification.size) + if not total: + return False + ground = int(np.count_nonzero(classification == GROUND_CLASS_CODE)) + return ground / total >= minimum_ratio diff --git a/B04_PreProcess/B04_PreProcess_Engine_Filter_Grid.py b/B04_PreProcess/B04_PreProcess_Engine_Filter_Grid.py index 3f0e0d0a..c96038de 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_Filter_Grid.py +++ b/B04_PreProcess/B04_PreProcess_Engine_Filter_Grid.py @@ -28,24 +28,21 @@ def filter_grid_min_z( if xyz.shape[0] == 0: return np.zeros(0, dtype=bool) - x_min, y_min, z_min = bounds[0, 0], bounds[1, 0], bounds[2, 0] + x_min, y_min = bounds[0, 0], bounds[1, 0] x_max, y_max = bounds[0, 1], bounds[1, 1] grid_width = int(np.ceil((x_max - x_min) / cell_size)) + 2 grid_height = int(np.ceil((y_max - y_min) / cell_size)) + 2 - minimum_z = np.full((grid_height, grid_width), np.inf, dtype=np.float32) + # float32로 담으면 셀 최저점(가장 확실한 지면점)이 반올림 때문에 자기 기준면보다 + # 낮아져 height_above < 0 으로 탈락한다 — 원본 좌표와 같은 float64로 둔다. + minimum_z = np.full((grid_height, grid_width), np.inf, dtype=np.float64) grid_x = np.clip(((xyz[:, 0] - x_min) / cell_size).astype(np.int64), 0, grid_width - 1) grid_y = np.clip(((xyz[:, 1] - y_min) / cell_size).astype(np.int64), 0, grid_height - 1) np.minimum.at(minimum_z, (grid_y, grid_x), xyz[:, 2]) - minimum_z[np.isinf(minimum_z)] = z_min - - try: - from scipy.ndimage import minimum_filter - - minimum_z = minimum_filter(minimum_z, size=3).astype(np.float32) - except ImportError: - pass + # 기준면은 점이 실제로 들어 있는 셀에서만 읽는다 — 3x3 최소필터로 이웃 셀 값을 + # 끌어오면 급경사에서 기준면이 경사만큼 파고들어(중앙 1.88m > 임계 1.5m) 지면점이 + # 통째로 탈락한다. 빈 셀은 아래 인덱싱에서 참조되지 않으므로 채울 필요가 없다. height_above = xyz[:, 2] - minimum_z[grid_y, grid_x] return np.asarray( (height_above >= 0.0) & (height_above <= height_threshold), diff --git a/B04_PreProcess/B04_PreProcess_Engine_Ground.py b/B04_PreProcess/B04_PreProcess_Engine_Ground.py index 041e0553..8500b8b4 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_Ground.py +++ b/B04_PreProcess/B04_PreProcess_Engine_Ground.py @@ -9,6 +9,10 @@ from typing import Any import numpy as np +from B04_PreProcess.B04_PreProcess_Engine_Filter_Classification import ( + filter_classification, + has_classified_ground, +) from B04_PreProcess.B04_PreProcess_Engine_Filter_CSF import filter_csf from B04_PreProcess.B04_PreProcess_Engine_Filter_Grid import filter_grid_min_z from B04_PreProcess.B04_PreProcess_Engine_Filter_PMF import filter_pmf @@ -20,13 +24,29 @@ _FILTERS = { "csf": filter_csf, "pmf": filter_pmf, "ransac": filter_ransac, + "classification": filter_classification, } +# LAS에 지면분류가 실제로 들어 있을 때만 쓰는 필터 (미분류 LAS에서는 마스크가 빈다) +_CONDITIONAL_FILTERS = {"classification": has_classified_ground} + def available_filters() -> tuple[str, ...]: return tuple(_FILTERS.keys()) +def detect_extra_filters( + structured_data: dict[str, Any] | np.lib.npyio.NpzFile, + filter_keys: tuple[str, ...] | list[str], +) -> list[str]: + """이 LAS에서만 쓸 수 있는 조건부 필터 중 아직 목록에 없는 것을 고른다.""" + return [ + key + for key, is_usable in _CONDITIONAL_FILTERS.items() + if key not in filter_keys and is_usable(structured_data) + ] + + def run_ground_filter( filter_key: str, structured_data: dict[str, Any] | np.lib.npyio.NpzFile ) -> np.ndarray: diff --git a/B04_PreProcess/B04_PreProcess_Schema.py b/B04_PreProcess/B04_PreProcess_Schema.py index cf6520ad..b13522af 100644 --- a/B04_PreProcess/B04_PreProcess_Schema.py +++ b/B04_PreProcess/B04_PreProcess_Schema.py @@ -9,7 +9,9 @@ from config.config_system import ( SURFACE_MODEL_SOURCE_FILTERS, ) -_ALLOWED_FILTERS = set(SURFACE_MODEL_SOURCE_FILTERS) | {"ransac"} +# classification은 LAS가 지면분류를 싣고 있을 때 엔진이 자동으로 붙이지만, +# 재계산 요청에서 이름으로 지정하는 것도 허용한다. +_ALLOWED_FILTERS = set(SURFACE_MODEL_SOURCE_FILTERS) | {"ransac", "classification"} _ALLOWED_METHODS = set(SURFACE_MODEL_PRECOMPUTE) diff --git a/config/config_system.py b/config/config_system.py index 9c310634..8bc08a47 100644 --- a/config/config_system.py +++ b/config/config_system.py @@ -106,7 +106,10 @@ SURFACE_CSF_CLOTH_RESOLUTION_M = float(os.getenv("SURFACE_CSF_CLOTH_RESOLUTION_M SURFACE_CSF_RIGIDNESS = int(os.getenv("SURFACE_CSF_RIGIDNESS", "1")) SURFACE_CSF_TIME_STEP = float(os.getenv("SURFACE_CSF_TIME_STEP", "0.65")) SURFACE_CSF_CLASS_THRESHOLD_M = float(os.getenv("SURFACE_CSF_CLASS_THRESHOLD_M", "0.5")) +# 낙하 후 스프링 정착에 쓰는 여유 반복 수 (실제 낙하 구간은 현장 기복에서 계산한다) SURFACE_CSF_ITERATIONS = int(os.getenv("SURFACE_CSF_ITERATIONS", "150")) +# 낙하+정착 반복 상한 — 기복 1500m(0.3185m/회)까지 커버한다 +SURFACE_CSF_MAX_ITERATIONS = int(os.getenv("SURFACE_CSF_MAX_ITERATIONS", "5000")) SURFACE_CSF_SLOPE_SMOOTH = os.getenv("SURFACE_CSF_SLOPE_SMOOTH", "True").lower() == "true" SURFACE_CSF_SLOPE_SMOOTH_THRESHOLD_M = float( os.getenv("SURFACE_CSF_SLOPE_SMOOTH_THRESHOLD_M", "1.8") @@ -126,6 +129,14 @@ SURFACE_RANSAC_ITERATIONS = int(os.getenv("SURFACE_RANSAC_ITERATIONS", "100")) SURFACE_RANSAC_LOCAL_GRID_SIZE_M = float(os.getenv("SURFACE_RANSAC_LOCAL_GRID_SIZE_M", "10.0")) SURFACE_RANSAC_SEED = int(os.getenv("SURFACE_RANSAC_SEED", "42")) +# LAS 자체 지면분류(class 2) 필터 — 이 비율 이상 분류돼 있을 때만 필터로 제공한다. +# 미분류 LAS는 classification이 전부 0이라 마스크가 비어 버린다. +SURFACE_CLASSIFIED_GROUND_MIN_RATIO = float( + os.getenv("SURFACE_CLASSIFIED_GROUND_MIN_RATIO", "0.002") +) +# 지면점 비율이 이 값 미만이면 필터가 사실상 실패한 것으로 보고 WARNING을 남긴다. +SURFACE_GROUND_RATIO_WARN = float(os.getenv("SURFACE_GROUND_RATIO_WARN", "0.01")) + # ───────────────────────────────────────────────────────────────────────── # 5-2. 지표면 모델 생성 파라미터 (TIN/DTM/NURBS/implicit/meshfree) # ─────────────────────────────────────────────────────────────────────────