Merge remote-tracking branch 'origin/sub_laptop_5' into main_laptop_1

This commit is contained in:
2026-11-01 23:03:44 +09:00
7 changed files with 318 additions and 11 deletions
@@ -0,0 +1,247 @@
"""B04 점진 TIN(PTD · Axelsson 꼴) 지면 필터 — 분류 없는 LAS 의 기본 필터.
1. 작은 칸(후보 칸) 최저점 = 후보 · 시작 간격 칸 최저점 = 씨앗.
2. 씨앗 삼각망 → 후보마다 삼각면까지 연직 거리 · 세 꼭짓점 각을 보고 문턱 안이면
삼각형마다 가장 가까운 하나를 지면에 더함 → 더 안 들어올 때까지 되풀이.
삼각망 밖 후보 · 격자는 맞닿은 테두리 삼각형 평면을 늘여 봄(밖으로 차례로 넓어짐).
굵게 → 가늘게: 후보가 모두 면 위로 뜬 칸(면이 처짐)에 반 간격 씨앗을 더해 다시 다짐.
끝에 둘레보다 최대 거리 넘게 튀거나 꺼진 외딴 꼭짓점(가시 · 잡점)을 걷음.
3. 마지막 지면 TIN 을 후보 칸 격자로 굳히고 모든 점을 면까지 연직 거리로 판정 —
면 아래 최대 거리 · 면 위 0.2 m(꼭짓점 = 칸 최저라 땅은 면 곁 · 위는 덤불).
장면 셋(DJI Terra 꼴) = 평지 3° · 0.3 m / 완경사 6° · 0.5 m / 급경사 10° · 0.7 m
· 기본 급경사(임도).
큰 파일 = 점을 덩이로 나눠 칸 최저만 모음 · 후보가 상한을 넘으면 후보 칸을 키움(메모리 상한).
"""
import math
from typing import Any
import numpy as np
from scipy import ndimage
from scipy.spatial import Delaunay, cKDTree
from config.config_system import SURFACE_PTD_SCENE, SURFACE_PTD_SEED_SPACING_M
# 장면 → (최대 각 °, 최대 거리 m)
PTD_SCENES = {"flat": (3.0, 0.3), "gentle": (6.0, 0.5), "steep": (10.0, 0.7)}
_CHUNK_POINTS = 5_000_000
def _cell_minima(xyz: np.ndarray, origin: tuple[float, float], cell: float, cols: int):
"""칸마다 가장 낮은 점의 (칸 번호, 점 번호)를 덩이로 나눠 모은다."""
keys, picks = [], []
for start in range(0, len(xyz), _CHUNK_POINTS):
part = xyz[start : start + _CHUNK_POINTS]
gx = ((part[:, 0] - origin[0]) / cell).astype(np.int64)
gy = ((part[:, 1] - origin[1]) / cell).astype(np.int64)
key = gy * cols + gx
order = np.lexsort((part[:, 2], key))
first = np.r_[True, key[order][1:] != key[order][:-1]]
keys.append(key[order][first])
picks.append(order[first] + start)
key = np.concatenate(keys)
pick = np.concatenate(picks)
order = np.lexsort((xyz[pick, 2], key))
first = np.r_[True, key[order][1:] != key[order][:-1]]
return key[order][first], pick[order][first]
def _lone_vertices(points: np.ndarray, gap: float, rounds: int = 10) -> np.ndarray:
"""둘레 꼭짓점 모두보다 gap 넘게 높거나(덤불 가시) 낮은(땅 밑 잡점) 꼭짓점 번호.
능선 · 골은 축 따라 비슷한 이웃이 있어 안 걸림."""
keep = np.arange(len(points))
for _ in range(rounds):
if len(keep) < 4:
break
indptr, nb = Delaunay(points[keep, :2]).vertex_neighbor_vertices
z = points[keep, 2]
low = np.minimum.reduceat(z[nb], indptr[:-1])
high = np.maximum.reduceat(z[nb], indptr[:-1])
lone = (z > high + gap) | (z < low - gap)
if not lone.any():
break
keep = keep[~lone]
return np.setdiff1d(np.arange(len(points)), keep)
def _lowest_per_cell(points: np.ndarray, origin: tuple[float, float], size: float):
"""칸마다 가장 낮은 점 — (칸 열, 칸 행, 점 번호)."""
gx = ((points[:, 0] - origin[0]) / size).astype(np.int64)
gy = ((points[:, 1] - origin[1]) / size).astype(np.int64)
key = gy * (gx.max() + 1) + gx
order = np.lexsort((points[:, 2], key))
first = np.r_[True, key[order][1:] != key[order][:-1]]
pick = order[first]
return gx[pick], gy[pick], pick
def _barycentric(tri: Delaunay, simplex: np.ndarray, q: np.ndarray) -> np.ndarray:
"""세 꼭짓점 무게(차례 = simplices 차례). tri.transform 은 BLAS 스레드 다툼으로
이 PC 에서 1 만 삼각형에 14 초 → 직접 셈."""
p = tri.points[tri.simplices[simplex]]
v0, v1, v2 = p[:, 1] - p[:, 0], p[:, 2] - p[:, 0], q - p[:, 0]
den = v0[:, 0] * v1[:, 1] - v1[:, 0] * v0[:, 1]
den = np.where(np.abs(den) < 1e-12, 1e-12, den)
l1 = (v2[:, 0] * v1[:, 1] - v1[:, 0] * v2[:, 1]) / den
l2 = (v0[:, 0] * v2[:, 1] - v2[:, 0] * v0[:, 1]) / den
return np.column_stack([1.0 - l1 - l2, l1, l2])
def _locate(tri: Delaunay, q: np.ndarray, max_steps: int = 200) -> np.ndarray:
"""점이 든 삼각형 번호 · 삼각망 밖이면 맞닿은 테두리 삼각형(그 평면으로 늘여 봄).
가장 가까운 꼭짓점에서 이웃으로 걸음 — scipy find_simplex 는 점마다 0 번에서 걸어
백만 점에 수십 분이 걸린다."""
_, nearest = cKDTree(tri.points).query(q)
simplex = tri.vertex_to_simplex[nearest].astype(np.int64)
active = np.arange(len(q))
for _ in range(max_steps):
if not len(active):
break
bary = _barycentric(tri, simplex[active], q[active])
worst = bary.argmin(axis=1)
moving = bary[np.arange(len(active)), worst] < -1e-9
step = tri.neighbors[simplex[active], worst]
walk = moving & (step >= 0)
simplex[active[walk]] = step[walk]
active = active[walk]
return simplex
def _plane_distance(tri: Delaunay, gz: np.ndarray, simplex: np.ndarray, q: np.ndarray):
"""점 q 의 삼각면까지 연직 거리와 세 꼭짓점에서 본 최대 각(°).
연직 = 비탈에서 수직 거리보다 관목을 덜 받음(수직은 45° 비탈 0.7 m 가 연직 1 m)."""
vertex = tri.simplices[simplex]
weight = _barycentric(tri, simplex, q[:, :2])
distance = np.abs(q[:, 2] - (weight * gz[vertex]).sum(axis=1))
reach = np.linalg.norm(q[:, None, :2] - tri.points[vertex], axis=2).min(axis=1)
angle = np.degrees(np.arctan(distance / np.maximum(reach, 1e-9)))
return distance, angle
def _tin_dz(ground: np.ndarray, q: np.ndarray) -> np.ndarray:
"""q 높이 − 지면 TIN 높이(연직)."""
tri = Delaunay(ground[:, :2])
simplex = _locate(tri, q[:, :2])
weight = _barycentric(tri, simplex, q[:, :2])
return q[:, 2] - (weight * ground[tri.simplices[simplex], 2]).sum(axis=1)
def _densify(
cand: np.ndarray, is_ground: np.ndarray, max_angle: float, max_distance: float, rounds: int
) -> None:
"""삼각형마다 문턱 안 가장 가까운 후보 하나를 지면에 더함 — 더 안 들어올 때까지."""
for _ in range(rounds):
ground = cand[is_ground]
tri = Delaunay(ground[:, :2])
pending = np.flatnonzero(~is_ground)
if not len(pending):
return
simplex = _locate(tri, cand[pending, :2])
distance, angle = _plane_distance(tri, ground[:, 2], simplex, cand[pending])
# 꼭짓점 곁 점은 각이 커지기만 함 → 최대 거리 절반 안이면 각을 안 봄
ok = (distance <= max_distance) & ((angle <= max_angle) | (distance <= max_distance / 2))
pending, simplex, distance = pending[ok], simplex[ok], distance[ok]
if not len(pending):
return
order = np.lexsort((distance, simplex))
first = np.r_[True, simplex[order][1:] != simplex[order][:-1]]
is_ground[pending[order[first]]] = True
if first.sum() < max(10, 0.001 * is_ground.sum()):
return
def filter_ptd(
structured_data: dict[str, Any] | np.lib.npyio.NpzFile,
scene: str = SURFACE_PTD_SCENE,
seed_spacing: float = SURFACE_PTD_SEED_SPACING_M,
max_angle: float | None = None,
max_distance: float | None = None,
above_distance: float = 0.2,
candidate_cell: float = 0.5,
max_candidates: int = 3_000_000,
max_iterations: int = 60,
refine_min_spacing: float = 2.5,
) -> np.ndarray:
"""점진 TIN 으로 지면 포인트의 불리언 마스크를 반환한다."""
if scene not in PTD_SCENES:
raise ValueError(f"점진 TIN 장면은 {', '.join(PTD_SCENES)} 중 하나여야 합니다.")
scene_angle, scene_distance = PTD_SCENES[scene]
max_angle = scene_angle if max_angle is None else float(max_angle)
max_distance = scene_distance if max_distance is None else float(max_distance)
for name, value in (("시작 간격", seed_spacing), ("후보 칸", candidate_cell)):
if not math.isfinite(value) or value <= 0:
raise ValueError(f"{name}은 0보다 큰 유한한 값이어야 합니다.")
if not (0 < max_angle < 90) or max_distance <= 0 or above_distance < 0:
raise ValueError("최대 각은 0~90°, 최대 거리는 0보다 커야 합니다.")
xyz = np.asarray(structured_data["xyz"], dtype=np.float64)
if xyz.ndim != 2 or xyz.shape[1] != 3:
raise ValueError("xyz 배열은 (N, 3) 형태여야 합니다.")
if len(xyz) < 3:
return np.ones(len(xyz), dtype=bool)
x_min, y_min = float(xyz[:, 0].min()), float(xyz[:, 1].min())
x_max, y_max = float(xyz[:, 0].max()), float(xyz[:, 1].max())
width, height = max(x_max - x_min, 1e-6), max(y_max - y_min, 1e-6)
# ponytail: 후보 상한은 칸을 키워 지킴(공간 타일 아님) · 수 km² 넘으면 겹침 타일로
cell = max(candidate_cell, math.sqrt(width * height / max_candidates))
cols = int(width / cell) + 1
cand_key, cand_idx = _cell_minima(xyz, (x_min, y_min), cell, cols)
cand = xyz[cand_idx]
# 씨앗 = 시작 간격 칸마다 가장 낮은 후보 + 자료 가장자리는 1/4 간격 칸마다(삼각망이
# 가장자리까지 닿게 — 밖은 늘여 볼 테두리 삼각형이 가늘어 땅을 놓침)
origin = (x_min, y_min)
is_ground = np.zeros(len(cand), dtype=bool)
is_ground[_lowest_per_cell(cand, origin, seed_spacing)[2]] = True
gx, gy, pick = _lowest_per_cell(cand, origin, seed_spacing / 4)
occupied = np.zeros((gy.max() + 1, gx.max() + 1), dtype=bool)
occupied[gy, gx] = True
edge = occupied & ~ndimage.binary_erosion(occupied, iterations=1, border_value=0)
is_ground[pick[edge[gy, gx]]] = True
_densify(cand, is_ground, max_angle, max_distance, max_iterations)
# 굵게 → 가늘게: 칸 안 후보가 모두 면보다 최대 거리 넘게 위면(땅이 모든 반사점 밑일 수
# 없음 = 성긴 씨앗 사이에서 면이 처짐) 면에 가장 가까운 후보를 씨앗으로 더하고 다시 다짐
spacing = seed_spacing / 2
while spacing >= refine_min_spacing:
dz = np.zeros(len(cand))
pending = np.flatnonzero(~is_ground)
dz[pending] = _tin_dz(cand[is_ground], cand[pending])
gx, gy, pick = _lowest_per_cell(np.column_stack([cand[:, :2], dz]), origin, spacing)
# 간격 절반보다 높이 뜬 칸 = 땅이 안 보이는 수관 · 덤불 → 씨앗으로 안 씀
sag = pick[(dz[pick] > max_distance) & (dz[pick] <= spacing / 2)]
if len(sag):
is_ground[sag] = True
_densify(cand, is_ground, max_angle, max_distance, max_iterations)
spacing /= 2
ground_idx = np.flatnonzero(is_ground)
is_ground[ground_idx[_lone_vertices(cand[ground_idx], max_distance)]] = False
# 마지막 지면 TIN → 후보 칸 격자 → 모든 점 연직 거리 판정
ground = cand[is_ground]
tri = Delaunay(ground[:, :2])
rows = int(height / cell) + 2
mesh_x, mesh_y = np.meshgrid(x_min + np.arange(cols + 1) * cell, y_min + np.arange(rows) * cell)
node = np.column_stack([mesh_x.ravel(), mesh_y.ravel()])
simplex = _locate(tri, node)
weight = _barycentric(tri, simplex, node)
grid = (weight * ground[tri.simplices[simplex], 2]).sum(axis=1).reshape(mesh_x.shape)
mask = np.zeros(len(xyz), dtype=bool)
for start in range(0, len(xyz), _CHUNK_POINTS):
part = xyz[start : start + _CHUNK_POINTS]
fx = np.clip((part[:, 0] - x_min) / cell, 0, cols - 1e-6)
fy = np.clip((part[:, 1] - y_min) / cell, 0, rows - 1 - 1e-6)
ix, iy = fx.astype(np.int64), fy.astype(np.int64)
tx, ty = fx - ix, fy - iy
z_tin = (
grid[iy, ix] * (1 - tx) * (1 - ty)
+ grid[iy, ix + 1] * tx * (1 - ty)
+ grid[iy + 1, ix] * (1 - tx) * ty
+ grid[iy + 1, ix + 1] * tx * ty
)
dz = part[:, 2] - z_tin
mask[start : start + len(part)] = (dz >= -max_distance) & (dz <= above_distance)
return mask
@@ -18,6 +18,7 @@ from B04_PreProcess.B04_PreProcess_Engine_Filter_Classification import (
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
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
@@ -27,6 +28,7 @@ _FILTERS = {
"csf": filter_csf,
"pmf": filter_pmf,
"ransac": filter_ransac,
"ptd": filter_ptd,
"classification": filter_classification,
}
+2 -2
View File
@@ -51,7 +51,7 @@ import { currentProjectId, guardProjectLeave } from "../A00_Common/b_project_swi
// 고를 수 있는 필터. 자동 전처리는 이 중 기본 하나만 만들고, 나머지는 관리자가
// 드롭다운에서 고를 때 그 조합만 계산한다 (2026-09-01 사용자 확정).
const SOURCE_FILTERS = ["classification", "grid_min_z", "csf", "pmf"] as const;
const SOURCE_FILTERS = ["classification", "ptd", "grid_min_z", "csf", "pmf"] as const;
const MODEL_METHODS = ["tin", "dtm", "nurbs", "implicit", "meshfree"] as const;
const DEFAULT_FILTER = "csf";
const DEFAULT_METHOD = "implicit"; // 확정 기본 표현(요구 558) — 확정본이 있으면 그것이 먼저
@@ -77,7 +77,7 @@ function buildSelectField(
label,
options: values.map((value) => ({
value,
text: value.replaceAll("_", " ").toUpperCase(),
text: value === "ptd" ? "점진 TIN" : value.replaceAll("_", " ").toUpperCase(),
})),
value: defaultValue,
size: 36,
@@ -31,6 +31,7 @@ SurfaceFilter = Literal[
"grid_min_z",
"csf",
"pmf",
"ptd",
"ransac",
"sheet_tin_sheet",
"sheet_tin",
+7 -2
View File
@@ -44,6 +44,11 @@ SURFACE_CSF_SLOPE_SMOOTH_THRESHOLD_M = float(
os.getenv("SURFACE_CSF_SLOPE_SMOOTH_THRESHOLD_M", "1.8")
)
# 점진 TIN(PTD) 지면 분류 — 장면 flat · gentle · steep(각 · 거리 묶음) · 시작 간격(씨앗 칸)
# 시작 간격 = 숲 10 m(DJI 기본 20 m 는 급한 볼록 비탈을 못 따라감)
SURFACE_PTD_SCENE = os.getenv("SURFACE_PTD_SCENE", "steep")
SURFACE_PTD_SEED_SPACING_M = float(os.getenv("SURFACE_PTD_SEED_SPACING_M", "10.0"))
# PMF (Progressive Morphological Filter) 지면 분류 파라미터
SURFACE_PMF_CELL_SIZE_M = float(os.getenv("SURFACE_PMF_CELL_SIZE_M", "2.0"))
SURFACE_PMF_MAX_WINDOW_SIZE = int(os.getenv("SURFACE_PMF_MAX_WINDOW_SIZE", "40"))
@@ -83,12 +88,12 @@ SURFACE_ROUTE_EDGE_TRIM_M = float(os.getenv("SURFACE_ROUTE_EDGE_TRIM_M", "3.0"))
# 고를 수 있는 지면 필터 전체 — 관리자가 B04 드롭다운에서 요청할 때 쓰인다.
# 자동 전처리가 미리 만드는 조합이 아니다(그건 SURFACE_AUTO_* 가 정한다).
SURFACE_MODEL_SOURCE_FILTERS = tuple(
os.getenv("SURFACE_MODEL_SOURCE_FILTERS", "classification,grid_min_z,csf,pmf").split(",")
os.getenv("SURFACE_MODEL_SOURCE_FILTERS", "classification,ptd,grid_min_z,csf,pmf").split(",")
)
# 자동 전처리가 만드는 조합 — 기본 필터 1종 × 아래 표현. 나머지는 관리자 요청 시 만든다.
SURFACE_AUTO_METHODS = tuple(os.getenv("SURFACE_AUTO_METHODS", "dtm").split(","))
# 기본 필터는 입력 LAS를 보고 정한다 — 지면분류(class 2)가 있으면 그것을, 없으면 아래 값.
SURFACE_AUTO_FALLBACK_FILTER = os.getenv("SURFACE_AUTO_FALLBACK_FILTER", "csf")
SURFACE_AUTO_FALLBACK_FILTER = os.getenv("SURFACE_AUTO_FALLBACK_FILTER", "ptd")
# dtm을 먼저 빌드해야 TIN 등고선 사전 캐시가 dtm footprint를 참조할 수 있다 (PLAN D-6)
SURFACE_MODEL_PRECOMPUTE = tuple(
os.getenv("SURFACE_MODEL_PRECOMPUTE", "dtm,tin,nurbs,implicit,meshfree").split(",")
+52
View File
@@ -0,0 +1,52 @@
"""642-2 점진 TIN(PTD) 지면 필터 지킴 — 합성 비탈 + 나무 점에서 땅 점 고름 · 급경사 놓침 상한."""
import numpy as np
import pytest
from B04_PreProcess.B04_PreProcess_Engine_Filter_PTD import filter_ptd
from B04_PreProcess.B04_PreProcess_Engine_Ground import available_filters, run_ground_filter
def _terrain(x, y):
# 서쪽 완경사(약 10 %) → 동쪽 급경사(끝 약 100 %) + 골 하나
return 0.1 * x + 0.0075 * np.maximum(x - 60, 0) ** 2 + 2.0 * np.sin(y / 15)
def _scene(seed=642):
rng = np.random.default_rng(seed)
n_ground = 60_000
gx, gy = rng.uniform(0, 120, n_ground), rng.uniform(0, 120, n_ground)
ground = np.column_stack([gx, gy, _terrain(gx, gy) + rng.normal(0, 0.03, n_ground)])
# 나무 = 줄기 둘레 수관(2~15 m 위) · 관목(0.9~2 m 위)
trunks = rng.uniform(0, 120, (120, 2))
tx = np.repeat(trunks[:, 0], 300) + rng.normal(0, 2.5, 36_000)
ty = np.repeat(trunks[:, 1], 300) + rng.normal(0, 2.5, 36_000)
tz = _terrain(tx, ty) + rng.uniform(2, 15, 36_000)
sx, sy = rng.uniform(0, 120, 8_000), rng.uniform(0, 120, 8_000)
sz = _terrain(sx, sy) + rng.uniform(0.9, 2.0, 8_000)
other = np.column_stack([np.r_[tx, sx], np.r_[ty, sy], np.r_[tz, sz]])
xyz = np.vstack([ground, other])
truth = np.r_[np.ones(n_ground, bool), np.zeros(len(other), bool)]
return xyz, truth
def test_ptd_picks_ground_on_slope_with_trees():
xyz, truth = _scene()
mask = filter_ptd({"xyz": xyz})
hit = (mask & truth).sum() / truth.sum()
wrong = (mask & ~truth).sum() / truth.sum()
assert hit >= 0.95, hit
assert wrong <= 0.02, wrong
# 급경사 칸(x > 90 · 기울기 55 % 넘음) 놓침 상한
steep = truth & (xyz[:, 0] > 90)
assert (~mask[steep]).mean() <= 0.08, (~mask[steep]).mean()
def test_ptd_registered_and_scenes():
assert "ptd" in available_filters()
xyz, _ = _scene(1)
flat = run_ground_filter("ptd", {"xyz": xyz})
assert flat.dtype == bool and len(flat) == len(xyz)
with pytest.raises(ValueError):
filter_ptd({"xyz": xyz}, scene="moon")
assert filter_ptd({"xyz": np.zeros((0, 3))}).shape == (0,)
+7 -7
View File
@@ -102,20 +102,20 @@ def test_classification_filter_reads_class_two() -> None:
assert mask.sum() == len(points[::4])
def test_auto_filters_prefer_classification_but_always_build_csf() -> None:
"""자동 전처리 목록 — 첫 항목이 기본 확정값. csf는 분류가 있어도 늘 함께 만든다."""
def test_auto_filters_prefer_classification_but_always_build_ptd() -> None:
"""자동 전처리 목록 — 첫 항목이 기본 확정값. ptd(점진 TIN)는 분류가 있어도 늘 함께 만든다."""
points = _sloped_ground(10.0)
unclassified = np.zeros(len(points), dtype=np.uint8)
classified = np.ones(len(points), dtype=np.uint8)
classified[::4] = 2
# 분류가 있으면 classification 이 기본이고, csf 는 검증용으로 같이 만든다.
# 분류가 있으면 classification 이 기본이고, ptd 는 검증용으로 같이 만든다.
assert resolve_auto_source_filters(_structured(points, classified)) == [
"classification",
"csf",
"ptd",
]
# 분류가 없으면 csf 하나뿐 — 빈 마스크로 모델을 만들지 않는다.
assert resolve_auto_source_filters(_structured(points, unclassified)) == ["csf"]
# 분류가 없으면 ptd 하나뿐 — 빈 마스크로 모델을 만들지 않는다.
assert resolve_auto_source_filters(_structured(points, unclassified)) == ["ptd"]
assert not has_classified_ground(_structured(points, unclassified))
# classification 이 없는 구조화 데이터(옛 캐시)에서도 죽지 않는다.
assert resolve_auto_source_filters(_structured(points)) == ["csf"]
assert resolve_auto_source_filters(_structured(points)) == ["ptd"]