398 lines
17 KiB
Python
398 lines
17 KiB
Python
from __future__ import annotations
|
|
|
|
import json
|
|
from pathlib import Path
|
|
from typing import Any
|
|
import numpy as np
|
|
from scipy import ndimage
|
|
from scipy.interpolate import RectBivariateSpline, RBFInterpolator
|
|
from skimage import measure
|
|
|
|
# config 파라미터 로드
|
|
import config
|
|
|
|
# 등고선 캐시 형식/추출 규칙이 바뀔 때 증가시킨다.
|
|
CONTOUR_EXTRACTOR_VERSION = 3
|
|
|
|
def extract_contours_from_grid(
|
|
x_coords: np.ndarray,
|
|
y_coords: np.ndarray,
|
|
z_grid: np.ndarray,
|
|
valid_mask: np.ndarray,
|
|
interval: float,
|
|
min_interval: float = 0.5,
|
|
scene_center: tuple[float, float, float] | None = None
|
|
) -> list[dict[str, Any]]:
|
|
"""
|
|
정규 표고 격자로부터 등고선 라인을 추출합니다.
|
|
scene_center가 (cx, cy, cz) 형태로 제공되면, 3D Scene 좌표계(Three.js)로 자동 변환하여 정점들을 저장합니다.
|
|
"""
|
|
interval = max(interval, min_interval)
|
|
|
|
# z_grid의 유효값 범위 획득
|
|
finite_mask = np.isfinite(z_grid)
|
|
if valid_mask is not None:
|
|
finite_mask &= valid_mask
|
|
|
|
if not finite_mask.any():
|
|
return []
|
|
|
|
z_min = float(np.min(z_grid[finite_mask]))
|
|
z_max = float(np.max(z_grid[finite_mask]))
|
|
|
|
# 등고선 생성 기준 레벨 결정
|
|
start_level = np.ceil(z_min / interval) * interval
|
|
levels = np.arange(start_level, z_max, interval)
|
|
if len(levels) == 0:
|
|
return []
|
|
|
|
# 최대 레벨 제한 (과도한 루프 방지)
|
|
if len(levels) > 500:
|
|
new_interval = (z_max - z_min) / 100.0
|
|
levels = np.arange(np.ceil(z_min / new_interval) * new_interval, z_max, new_interval)
|
|
interval = new_interval
|
|
|
|
contours_geojson_list = []
|
|
|
|
# marching squares의 NaN 문제를 예방하기 위해,
|
|
# 마스크가 아닌 영역(~valid_mask)을 z_min보다 충분히 낮은 sentinel(z_min - 1000)으로 채웁니다.
|
|
z_grid_masked = z_grid.copy()
|
|
if valid_mask is not None:
|
|
z_grid_masked[~valid_mask] = z_min - 1000.0
|
|
# 유효한 영역 내의 NaN 값들도 함께 채워줍니다.
|
|
invalid_mask = ~np.isfinite(z_grid_masked)
|
|
if invalid_mask.any():
|
|
z_grid_masked[invalid_mask] = z_min - 1000.0
|
|
|
|
cx, cy, cz = scene_center if scene_center is not None else (0.0, 0.0, 0.0)
|
|
|
|
# marching squares로 레벨별 등고선 추출
|
|
for level in levels:
|
|
contours = measure.find_contours(z_grid_masked, level)
|
|
for contour in contours:
|
|
current_segment = []
|
|
for y_idx, x_idx in contour:
|
|
# 격자 경계 클램프
|
|
x_idx_clamped = np.clip(x_idx, 0, len(x_coords) - 1)
|
|
y_idx_clamped = np.clip(y_idx, 0, len(y_coords) - 1)
|
|
|
|
# 인접 정수 인덱스
|
|
x_0 = int(np.floor(x_idx_clamped))
|
|
x_1 = int(np.ceil(x_idx_clamped))
|
|
y_0 = int(np.floor(y_idx_clamped))
|
|
y_1 = int(np.ceil(y_idx_clamped))
|
|
|
|
# 최종 정점이 유효 마스크 내부인지 엄격히 판별
|
|
is_valid = True
|
|
if valid_mask is not None:
|
|
if not (valid_mask[y_0, x_0] and valid_mask[y_0, x_1] and valid_mask[y_1, x_0] and valid_mask[y_1, x_1]):
|
|
is_valid = False
|
|
|
|
if not is_valid:
|
|
# 유효하지 않은 영역을 만났다면, 현재까지 세그먼트가 존재할 경우 저장하고 새 세그먼트 준비
|
|
if len(current_segment) >= 2:
|
|
mid_idx = len(current_segment) // 2
|
|
contours_geojson_list.append({
|
|
"level": float(level),
|
|
"coordinates": current_segment,
|
|
"label_position": current_segment[mid_idx]
|
|
})
|
|
current_segment = []
|
|
continue
|
|
|
|
# 소수부 비율
|
|
tx = x_idx_clamped - x_0
|
|
ty = y_idx_clamped - y_0
|
|
|
|
# 보간한 X, Y 좌표 계산
|
|
x_val = (1.0 - tx) * x_coords[x_0] + tx * x_coords[x_1]
|
|
y_val = (1.0 - ty) * y_coords[y_0] + ty * y_coords[y_1]
|
|
|
|
if scene_center is not None:
|
|
# Three.js 씬 좌표 매핑 변환 적용
|
|
pt_x = round(float(x_val - cx), 3)
|
|
pt_y = round(float(level - cz), 3)
|
|
pt_z = round(float(-(y_val - cy)), 3)
|
|
current_segment.append([pt_x, pt_y, pt_z])
|
|
else:
|
|
current_segment.append([round(float(x_val), 3), round(float(y_val), 3), round(float(level), 3)])
|
|
|
|
# 마지막 남은 세그먼트가 유효하다면 추가
|
|
if len(current_segment) >= 2:
|
|
mid_idx = len(current_segment) // 2
|
|
contours_geojson_list.append({
|
|
"level": float(level),
|
|
"coordinates": current_segment,
|
|
"label_position": current_segment[mid_idx]
|
|
})
|
|
|
|
return contours_geojson_list
|
|
|
|
def _load_footprint_mask(
|
|
model_npz_path: Path,
|
|
x_coords: np.ndarray,
|
|
y_coords: np.ndarray,
|
|
) -> np.ndarray | None:
|
|
"""같은 source filter의 DTM valid_mask(실제 지면 footprint)를 현재 method 격자에
|
|
최근접 리샘플하여 반환한다. DTM 캐시가 없으면 None.
|
|
(TIN/meshfree의 convex-hull 누출, NURBS/implicit의 직사각형 영역 누출을 모두 차단)
|
|
"""
|
|
stem = Path(model_npz_path).stem # 예: "tin_grid_min_z_smooth" 또는 "tin_grid_min_z"
|
|
if stem.endswith("_smooth"):
|
|
stem = stem[:-7]
|
|
parts = stem.split("_", 1)
|
|
if len(parts) < 2:
|
|
return None
|
|
filter_key = parts[1] # "grid_min_z", "csf", "pmf"
|
|
dtm_path = Path(model_npz_path).parent / f"dtm_{filter_key}.npz"
|
|
if not dtm_path.exists():
|
|
return None
|
|
try:
|
|
d = np.load(dtm_path)
|
|
dtm_x = np.asarray(d["x"]).ravel()
|
|
dtm_y = np.asarray(d["y"]).ravel()
|
|
dtm_mask = np.asarray(d["valid_mask"], dtype=bool)
|
|
except Exception:
|
|
return None
|
|
if len(dtm_x) < 2 or len(dtm_y) < 2:
|
|
return None
|
|
|
|
def _nearest_idx(axis: np.ndarray, coords: np.ndarray) -> np.ndarray:
|
|
ascending = bool(axis[0] <= axis[-1])
|
|
a = axis if ascending else axis[::-1]
|
|
idx = np.clip(np.searchsorted(a, coords), 1, len(a) - 1)
|
|
idx = np.where(np.abs(a[idx - 1] - coords) <= np.abs(a[idx] - coords), idx - 1, idx)
|
|
return idx if ascending else (len(axis) - 1 - idx)
|
|
|
|
xi = _nearest_idx(dtm_x, np.asarray(x_coords, dtype=np.float64))
|
|
yi = _nearest_idx(dtm_y, np.asarray(y_coords, dtype=np.float64))
|
|
return dtm_mask[np.ix_(yi, xi)]
|
|
|
|
def _apply_footprint(
|
|
model_npz_path: Path,
|
|
x_coords: np.ndarray,
|
|
y_coords: np.ndarray,
|
|
valid_mask: np.ndarray,
|
|
) -> np.ndarray:
|
|
"""현재 valid_mask에 DTM 기반 footprint 마스크를 교집합으로 적용한다. 형상이 약간 다르면 zoom/nearest로 자동 형상 조절하여 강제 적용한다."""
|
|
fp = _load_footprint_mask(model_npz_path, x_coords, y_coords)
|
|
if fp is not None:
|
|
if fp.shape == valid_mask.shape:
|
|
return valid_mask & fp
|
|
else:
|
|
# 해상도가 다를 경우 최근접 리샘플링하여 크기를 맞춤
|
|
from scipy.ndimage import zoom
|
|
zoom_y = valid_mask.shape[0] / fp.shape[0]
|
|
zoom_x = valid_mask.shape[1] / fp.shape[1]
|
|
fp_resized = zoom(fp.astype(float), (zoom_y, zoom_x), order=0) > 0.5
|
|
if fp_resized.shape == valid_mask.shape:
|
|
return valid_mask & fp_resized
|
|
return valid_mask
|
|
|
|
|
|
def _tin_face_coverage_mask(
|
|
vertices: np.ndarray,
|
|
faces: np.ndarray,
|
|
xx: np.ndarray,
|
|
yy: np.ndarray,
|
|
) -> np.ndarray:
|
|
"""저장된 TIN 면이 실제로 덮는 XY 영역만 True로 반환한다.
|
|
|
|
정점으로 Delaunay를 다시 수행하면 TIN 생성 단계에서 제거한 긴 면과
|
|
footprint 경계 면이 되살아날 수 있다. 따라서 실제 저장된 faces의 경계
|
|
에지를 polygonize하여 렌더링 메시와 동일한 coverage를 사용한다.
|
|
"""
|
|
vertices = np.asarray(vertices)
|
|
faces = np.asarray(faces, dtype=np.int64)
|
|
if vertices.ndim != 2 or vertices.shape[1] < 2 or not len(faces):
|
|
return np.zeros(xx.shape, dtype=bool)
|
|
|
|
# 삼각형마다 세 변을 모은 뒤 한 번만 나타나는 변을 메시 경계로 본다.
|
|
edges = np.vstack((faces[:, [0, 1]], faces[:, [1, 2]], faces[:, [2, 0]]))
|
|
edges = np.sort(edges, axis=1)
|
|
unique_edges, counts = np.unique(edges, axis=0, return_counts=True)
|
|
boundary_edges = unique_edges[counts == 1]
|
|
if not len(boundary_edges):
|
|
return np.zeros(xx.shape, dtype=bool)
|
|
|
|
from shapely import get_parts, intersects_xy, linestrings, polygonize, union_all
|
|
|
|
boundary_lines = linestrings(vertices[boundary_edges, :2])
|
|
polygons = list(get_parts(polygonize(boundary_lines)))
|
|
if not polygons:
|
|
return np.zeros(xx.shape, dtype=bool)
|
|
coverage = union_all(polygons)
|
|
|
|
# shapely.intersects_xy의 x, y 인자는 1차원 float64 배열이어야 하므로 평탄화 후 타입 캐스팅 적용
|
|
xx_flat = np.asarray(xx, dtype=np.float64).ravel()
|
|
yy_flat = np.asarray(yy, dtype=np.float64).ravel()
|
|
|
|
res_flat = np.asarray(intersects_xy(coverage, xx_flat, yy_flat), dtype=bool)
|
|
return res_flat.reshape(xx.shape)
|
|
|
|
def extract_contours(
|
|
model_npz_path: Path,
|
|
representation: str,
|
|
interval: float,
|
|
target_grid_m: float = 1.0,
|
|
scene_center: tuple[float, float, float] | None = None
|
|
) -> list[dict[str, Any]]:
|
|
"""
|
|
15종 표현형식의 npz 모델 파일로부터 표고 격자를 환원한 후,
|
|
정해진 interval 간격의 등고선 리스트를 추출합니다.
|
|
"""
|
|
model_npz_path = Path(model_npz_path)
|
|
if not model_npz_path.exists():
|
|
raise FileNotFoundError(f"모델 파일이 존재하지 않습니다: {model_npz_path}")
|
|
|
|
data = np.load(model_npz_path)
|
|
|
|
# 1. 표현형식별 정규 격자(elevation grid) 환원
|
|
if representation == "regular_grid":
|
|
# DTM 격자 다운샘플
|
|
x_coords = data["x"]
|
|
y_coords = data["y"]
|
|
z_grid = data["z"]
|
|
valid_mask = data["valid_mask"]
|
|
|
|
# dtm의 dtm_grid_resolution_meters는 보통 0.05m이므로 1m 등으로 다운샘플
|
|
current_res = (x_coords[-1] - x_coords[0]) / (len(x_coords) - 1)
|
|
step = max(1, int(round(target_grid_m / current_res)))
|
|
|
|
if step > 1:
|
|
x_coords_sub = x_coords[::step]
|
|
y_coords_sub = y_coords[::step]
|
|
z_grid_sub = z_grid[::step, ::step]
|
|
valid_mask_sub = valid_mask[::step, ::step]
|
|
return extract_contours_from_grid(x_coords_sub, y_coords_sub, z_grid_sub, valid_mask_sub, interval, scene_center=scene_center)
|
|
else:
|
|
return extract_contours_from_grid(x_coords, y_coords, z_grid, valid_mask, interval, scene_center=scene_center)
|
|
|
|
elif representation == "triangular_mesh":
|
|
# TIN 메쉬 보간
|
|
vertices = data["vertices"] # (N, 3)
|
|
faces = data["faces"]
|
|
|
|
# 1.0m 격자로 샘플링
|
|
x_min, x_max = float(np.min(vertices[:, 0])), float(np.max(vertices[:, 0]))
|
|
y_min, y_max = float(np.min(vertices[:, 1])), float(np.max(vertices[:, 1]))
|
|
|
|
cols = max(2, int(np.ceil((x_max - x_min) / target_grid_m)) + 1)
|
|
rows = max(2, int(np.ceil((y_max - y_min) / target_grid_m)) + 1)
|
|
|
|
x_coords = np.linspace(x_min, x_max, cols, dtype=np.float32)
|
|
y_coords = np.linspace(y_min, y_max, rows, dtype=np.float32)
|
|
xx, yy = np.meshgrid(x_coords, y_coords)
|
|
|
|
# scipy.interpolate.griddata를 사용해 TIN 지면 보간
|
|
from scipy.interpolate import griddata
|
|
z_grid = griddata(vertices[:, :2], vertices[:, 2], (xx, yy), method="linear")
|
|
|
|
# 중요: vertices를 다시 Delaunay하면 생성 단계에서 제거한 삼각형이
|
|
# 복원된다. 저장된 faces의 실제 XY coverage만 유효 영역으로 사용한다.
|
|
face_mask = _tin_face_coverage_mask(vertices, faces, xx, yy)
|
|
valid_mask = np.isfinite(z_grid) & face_mask
|
|
|
|
# 실제 지면 footprint와 교집합 (DTM 기준)
|
|
valid_mask = _apply_footprint(model_npz_path, x_coords, y_coords, valid_mask)
|
|
|
|
return extract_contours_from_grid(x_coords, y_coords, z_grid, valid_mask, interval, scene_center=scene_center)
|
|
|
|
elif representation == "bspline_surface":
|
|
# NURBS 제어점 로드 및 평가
|
|
control_x = data["control_x"]
|
|
control_y = data["control_y"]
|
|
control_z = data["control_z"]
|
|
degree = int(data["degree"][0])
|
|
|
|
spline = RectBivariateSpline(
|
|
control_y, control_x, control_z,
|
|
kx=min(degree, len(control_y) - 1),
|
|
ky=min(degree, len(control_x) - 1),
|
|
s=float(len(control_x) * len(control_y)) * 0.01
|
|
)
|
|
|
|
# 1.0m 격자로 평가
|
|
x_min, x_max = float(control_x[0]), float(control_x[-1])
|
|
y_min, y_max = float(control_y[0]), float(control_y[-1])
|
|
cols = max(2, int(np.ceil((x_max - x_min) / target_grid_m)) + 1)
|
|
rows = max(2, int(np.ceil((y_max - y_min) / target_grid_m)) + 1)
|
|
|
|
x_coords = np.linspace(x_min, x_max, cols, dtype=np.float32)
|
|
y_coords = np.linspace(y_min, y_max, rows, dtype=np.float32)
|
|
|
|
z_grid = np.asarray(spline(y_coords, x_coords), dtype=np.float32)
|
|
# NURBS는 직사각형 도메인 전체를 평가하므로 실제 지면 footprint로 제한한다.
|
|
valid_mask = np.ones_like(z_grid, dtype=bool)
|
|
valid_mask = _apply_footprint(model_npz_path, x_coords, y_coords, valid_mask)
|
|
|
|
return extract_contours_from_grid(x_coords, y_coords, z_grid, valid_mask, interval, scene_center=scene_center)
|
|
|
|
elif representation == "local_rbf_height_field":
|
|
# Implicit RBF 평가
|
|
centers_xy = data["centers_xy"]
|
|
center_z = data["center_z"]
|
|
smoothing = float(data["smoothing"][0])
|
|
|
|
interpolator = RBFInterpolator(
|
|
centers_xy.astype(np.float64),
|
|
center_z.astype(np.float64),
|
|
neighbors=min(64, len(centers_xy)),
|
|
smoothing=smoothing,
|
|
kernel="thin_plate_spline"
|
|
)
|
|
|
|
x_min, x_max = float(np.min(centers_xy[:, 0])), float(np.max(centers_xy[:, 0]))
|
|
y_min, y_max = float(np.min(centers_xy[:, 1])), float(np.max(centers_xy[:, 1]))
|
|
cols = max(2, int(np.ceil((x_max - x_min) / target_grid_m)) + 1)
|
|
rows = max(2, int(np.ceil((y_max - y_min) / target_grid_m)) + 1)
|
|
|
|
x_coords = np.linspace(x_min, x_max, cols, dtype=np.float32)
|
|
y_coords = np.linspace(y_min, y_max, rows, dtype=np.float32)
|
|
xx, yy = np.meshgrid(x_coords, y_coords)
|
|
|
|
query = np.column_stack([xx.ravel(), yy.ravel()])
|
|
z_values = interpolator(query).astype(np.float32)
|
|
z_grid = z_values.reshape(len(y_coords), len(x_coords))
|
|
# Implicit RBF도 직사각형 도메인 전체를 평가하므로 실제 지면 footprint로 제한한다.
|
|
valid_mask = np.ones_like(z_grid, dtype=bool)
|
|
valid_mask = _apply_footprint(model_npz_path, x_coords, y_coords, valid_mask)
|
|
|
|
return extract_contours_from_grid(x_coords, y_coords, z_grid, valid_mask, interval, scene_center=scene_center)
|
|
|
|
elif representation == "meshfree_surfels":
|
|
# Meshfree 포인트 데이터 보간
|
|
points = data["points"] # (N, 3)
|
|
|
|
x_min, x_max = float(np.min(points[:, 0])), float(np.max(points[:, 0]))
|
|
y_min, y_max = float(np.min(points[:, 1])), float(np.max(points[:, 1]))
|
|
cols = max(2, int(np.ceil((x_max - x_min) / target_grid_m)) + 1)
|
|
rows = max(2, int(np.ceil((y_max - y_min) / target_grid_m)) + 1)
|
|
|
|
x_coords = np.linspace(x_min, x_max, cols, dtype=np.float32)
|
|
y_coords = np.linspace(y_min, y_max, rows, dtype=np.float32)
|
|
xx, yy = np.meshgrid(x_coords, y_coords)
|
|
|
|
from scipy.interpolate import griddata
|
|
from scipy.spatial import Delaunay
|
|
z_grid = griddata(points[:, :2], points[:, 2], (xx, yy), method="linear")
|
|
|
|
valid_mask = np.isfinite(z_grid)
|
|
|
|
# Meshfree도 griddata(linear)라 TIN과 동일하게 convex hull 경계 마스킹
|
|
try:
|
|
tri = Delaunay(points[:, :2])
|
|
query_pts = np.column_stack([xx.ravel(), yy.ravel()])
|
|
hull_inside = tri.find_simplex(query_pts) >= 0
|
|
hull_mask = hull_inside.reshape(xx.shape)
|
|
valid_mask = valid_mask & hull_mask
|
|
except Exception:
|
|
pass
|
|
|
|
valid_mask = _apply_footprint(model_npz_path, x_coords, y_coords, valid_mask)
|
|
|
|
return extract_contours_from_grid(x_coords, y_coords, z_grid, valid_mask, interval, scene_center=scene_center)
|
|
|
|
else:
|
|
raise ValueError(f"지원하지 않는 표현 방식입니다: {representation}")
|