eomsangdonandClaude Opus 5
9a5496c47f
csf-nurbs 등고선 캐시가 Maximum allowed size exceeded로 실패하던 문제. 평활계수 s가 제어점당 잔차제곱 0.01(RMS 0.1m)로 너무 빡빡해 FITPACK이 제어점 수(122x113)보다 많은 knot(126x117)을 밀어넣고 계수가 발산했다. 실측 표고: csf 1e105, grid_min_z 1e10. grid_min_z는 footprint 안쪽만 보고 z 범위를 잡아 조용히 통과하고 있었다. - ModelContext에 fit_nurbs_spline/evaluate_nurbs_spline 신설. 모델 빌더와 등고선 엔진이 같은 곡면을 각각 만들면서 s 식이 두 곳에 중복돼 있었다. - NURBS_RESIDUAL_PER_CONTROL_POINT = 1.0 (RMS 1m). knot이 17~35로 떨어지고 세 지면 필터 모두 데이터 표고 범위 안에 머문다. NURBS는 평활 곡면 표현이라 이 정도 완화가 맞다. - 평가 결과가 제어 표고 범위(±50% 여유)를 벗어나면 잘라내고 경고를 남긴다. 방치하면 float32 캐스팅에서 inf가 되어 프리뷰 색상(nan)과 등고선까지 번진다. - extract_contours_from_grid: 레벨 수 상한 검사를 np.arange 앞으로 옮겼다. 뒤에 있어서 상한이 무용지물이었다 — 배열이 만들어지기 전에 터진다. 확인: 3필터 모두 곡면이 제어 표고 범위 내(grid_min_z 503.6~573.1, csf 504.7~552.3, pmf 504.6~551.8), 등고선 0.55~1.64초에 정상 산출. 표고 범위 1e12 격자도 예외 없이 처리. tmp/tests 59건 통과. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
472 lines
20 KiB
Python
472 lines
20 KiB
Python
"""B04 등고선 추출 엔진.
|
|
|
|
5종 표현(regular_grid/triangular_mesh/bspline_surface/local_rbf_height_field/
|
|
meshfree_surfels)의 npz 모델에서 표고 격자를 환원하고, marching squares로
|
|
지정 간격 등고선 라인을 추출한다. DTM valid_mask를 footprint로 사용해
|
|
경계 누출을 차단한다.
|
|
"""
|
|
|
|
from pathlib import Path
|
|
from typing import Any
|
|
|
|
import numpy as np
|
|
from scipy.interpolate import RBFInterpolator
|
|
from skimage import measure
|
|
|
|
# 등고선 캐시 형식/추출 규칙이 바뀔 때 증가시킨다.
|
|
CONTOUR_EXTRACTOR_VERSION = 4
|
|
|
|
|
|
def extract_contours_from_grid(
|
|
x_coords: np.ndarray,
|
|
y_coords: np.ndarray,
|
|
z_grid: np.ndarray,
|
|
valid_mask: np.ndarray | None,
|
|
interval: float,
|
|
min_interval: float = 0.5,
|
|
scene_center: tuple[float, float, float] | None = None,
|
|
) -> list[dict[str, Any]]:
|
|
"""정규 표고 격자로부터 등고선 라인을 추출한다."""
|
|
interval = max(interval, min_interval)
|
|
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]))
|
|
span = z_max - z_min
|
|
if span <= 0:
|
|
return []
|
|
# 레벨 수 상한은 arange **앞에서** 건다. 뒤에 두면 표고 범위가 비정상적으로 넓을 때
|
|
# np.arange가 배열을 만들다 "Maximum allowed size exceeded"로 먼저 터진다.
|
|
if span / interval > 500:
|
|
interval = span / 100.0
|
|
start_level = np.ceil(z_min / interval) * interval
|
|
levels = np.arange(start_level, z_max, interval)
|
|
if len(levels) == 0:
|
|
return []
|
|
|
|
contours_geojson_list: list[dict[str, Any]] = []
|
|
|
|
# marching squares의 NaN 문제 예방: 무효 영역을 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
|
|
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)
|
|
|
|
for level in levels:
|
|
for contour in measure.find_contours(z_grid_masked, level):
|
|
current_segment: list[list[float]] = []
|
|
for y_idx, x_idx in contour:
|
|
x_idx_c = np.clip(x_idx, 0, len(x_coords) - 1)
|
|
y_idx_c = np.clip(y_idx, 0, len(y_coords) - 1)
|
|
x_0, x_1 = int(np.floor(x_idx_c)), int(np.ceil(x_idx_c))
|
|
y_0, y_1 = int(np.floor(y_idx_c)), int(np.ceil(y_idx_c))
|
|
|
|
is_valid = True
|
|
if valid_mask is not None and 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_c - x_0
|
|
ty = y_idx_c - y_0
|
|
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:
|
|
current_segment.append(
|
|
[
|
|
round(float(x_val - cx), 3),
|
|
round(float(level - cz), 3),
|
|
round(float(-(y_val - cy)), 3),
|
|
]
|
|
)
|
|
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를 현재 격자에 최근접 리샘플한다."""
|
|
stem = Path(model_npz_path).stem
|
|
if stem.endswith("_smooth"):
|
|
stem = stem[:-7]
|
|
parts = stem.split("_", 1)
|
|
if len(parts) < 2:
|
|
return None
|
|
filter_key = parts[1]
|
|
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를 교집합으로 적용한다 (형상 다르면 최근접 리샘플)."""
|
|
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
|
|
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로 반환한다."""
|
|
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)
|
|
|
|
import affine
|
|
import rasterio.features
|
|
from shapely import get_parts, linestrings, polygonize
|
|
|
|
boundary_lines = linestrings(vertices[boundary_edges, :2])
|
|
polygons = list(get_parts(polygonize(boundary_lines)))
|
|
if not polygons:
|
|
return np.zeros(xx.shape, dtype=bool)
|
|
|
|
x_coords = xx[0, :]
|
|
y_coords = yy[:, 0]
|
|
dx = float(x_coords[1] - x_coords[0]) if len(x_coords) > 1 else 1.0
|
|
dy = float(y_coords[1] - y_coords[0]) if len(y_coords) > 1 else 1.0
|
|
|
|
transform = affine.Affine(dx, 0.0, x_coords[0] - dx / 2.0, 0.0, dy, y_coords[0] - dy / 2.0)
|
|
|
|
mask = rasterio.features.rasterize(
|
|
polygons, out_shape=xx.shape, transform=transform, fill=0, default_value=1, dtype="uint8"
|
|
)
|
|
return mask.astype(bool)
|
|
|
|
|
|
def _grid_axes(x_min: float, x_max: float, y_min: float, y_max: float, target_grid_m: float):
|
|
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)
|
|
return x_coords, y_coords
|
|
|
|
|
|
# 후보 셀 배열이 한 번에 커지지 않도록 삼각형을 나눠 처리하는 단위.
|
|
_RASTER_CHUNK = 200_000
|
|
# 격자점이 변 위에 놓였을 때 양쪽 삼각형 모두에서 탈락하지 않도록 두는 여유.
|
|
_BARYCENTRIC_EPS = 1e-9
|
|
|
|
|
|
def rasterize_triangle_mesh(
|
|
vertices: np.ndarray,
|
|
triangles: np.ndarray,
|
|
x_coords: np.ndarray,
|
|
y_coords: np.ndarray,
|
|
) -> np.ndarray:
|
|
"""삼각망을 규칙격자에 직접 래스터화해 표고 격자를 만든다.
|
|
|
|
삼각형마다 자기가 덮는 격자 셀만 무게중심 선형보간으로 채우므로, 질의점별
|
|
위치탐색(scipy `find_simplex`)이 사라진다. 삼각망 밖 셀은 NaN으로 남아
|
|
기존 `griddata` 결과와 같은 유효 판정을 유지한다.
|
|
"""
|
|
rows, cols = len(y_coords), len(x_coords)
|
|
z_grid = np.full((rows, cols), np.nan, dtype=np.float64)
|
|
vertices = np.asarray(vertices, dtype=np.float64)
|
|
triangles = np.asarray(triangles, dtype=np.int64)
|
|
if not len(triangles) or vertices.ndim != 2 or vertices.shape[1] < 3:
|
|
return z_grid
|
|
|
|
axis_x = np.asarray(x_coords, dtype=np.float64)
|
|
axis_y = np.asarray(y_coords, dtype=np.float64)
|
|
|
|
def _index_range(axis: np.ndarray, low: np.ndarray, high: np.ndarray, size: int):
|
|
"""축 배열을 직접 탐색해 [low, high]를 덮는 인덱스 구간을 낸다.
|
|
|
|
`(좌표 - 원점) / 간격` 식은 쓰지 않는다. `_grid_axes`가 float32 축을 만드는데
|
|
UTM 좌표(18만대)에서 float32 해상도가 0.015625m라 격자 간격이 0.984~1.0으로
|
|
흔들리고, 등간격으로 환산하면 최대 0.81셀까지 어긋난다. 삼각형이 평균 0.3m라
|
|
그만큼 어긋나면 격자점을 통째로 놓친다 (2026-08-17 실측).
|
|
"""
|
|
ascending = size < 2 or axis[-1] >= axis[0]
|
|
sorted_axis = axis if ascending else axis[::-1]
|
|
lo = np.searchsorted(sorted_axis, low, side="left")
|
|
hi = np.searchsorted(sorted_axis, high, side="right") - 1
|
|
if ascending:
|
|
return lo, hi
|
|
return size - 1 - hi, size - 1 - lo
|
|
|
|
for begin in range(0, len(triangles), _RASTER_CHUNK):
|
|
chunk = triangles[begin : begin + _RASTER_CHUNK]
|
|
corners = vertices[chunk]
|
|
ax, ay, az = corners[:, 0, 0], corners[:, 0, 1], corners[:, 0, 2]
|
|
bx, by, bz = corners[:, 1, 0], corners[:, 1, 1], corners[:, 1, 2]
|
|
cx, cy, cz = corners[:, 2, 0], corners[:, 2, 1], corners[:, 2, 2]
|
|
|
|
denom = (by - cy) * (ax - cx) + (cx - bx) * (ay - cy)
|
|
col_lo, col_hi = _index_range(
|
|
axis_x, np.minimum(np.minimum(ax, bx), cx), np.maximum(np.maximum(ax, bx), cx), cols
|
|
)
|
|
row_lo, row_hi = _index_range(
|
|
axis_y, np.minimum(np.minimum(ay, by), cy), np.maximum(np.maximum(ay, by), cy), rows
|
|
)
|
|
|
|
# 클립 전에 격자 밖 삼각형을 걸러야 경계에 눌러붙지 않는다.
|
|
selected = (
|
|
np.isfinite(denom)
|
|
& (np.abs(denom) > 1e-12)
|
|
& (col_hi >= col_lo)
|
|
& (row_hi >= row_lo)
|
|
& (col_hi >= 0)
|
|
& (row_hi >= 0)
|
|
& (col_lo <= cols - 1)
|
|
& (row_lo <= rows - 1)
|
|
)
|
|
if not selected.any():
|
|
continue
|
|
|
|
col_lo = np.clip(col_lo[selected], 0, cols - 1).astype(np.int64)
|
|
col_hi = np.clip(col_hi[selected], 0, cols - 1).astype(np.int64)
|
|
row_lo = np.clip(row_lo[selected], 0, rows - 1).astype(np.int64)
|
|
row_hi = np.clip(row_hi[selected], 0, rows - 1).astype(np.int64)
|
|
width = col_hi - col_lo + 1
|
|
counts = width * (row_hi - row_lo + 1)
|
|
total = int(counts.sum())
|
|
if total <= 0:
|
|
continue
|
|
|
|
# 삼각형별 후보 셀 수가 제각각이라 반복 인덱스를 펼쳐서 한 번에 계산한다.
|
|
starts = np.concatenate(([0], np.cumsum(counts)[:-1]))
|
|
local = np.arange(total, dtype=np.int64) - np.repeat(starts, counts)
|
|
width_rep = np.repeat(width, counts)
|
|
offset_row = local // width_rep
|
|
cell_row = np.repeat(row_lo, counts) + offset_row
|
|
cell_col = np.repeat(col_lo, counts) + (local - offset_row * width_rep)
|
|
|
|
point_x = axis_x[cell_col]
|
|
point_y = axis_y[cell_row]
|
|
ax_r, ay_r, az_r = (
|
|
np.repeat(ax[selected], counts),
|
|
np.repeat(ay[selected], counts),
|
|
np.repeat(az[selected], counts),
|
|
)
|
|
bx_r, by_r, bz_r = (
|
|
np.repeat(bx[selected], counts),
|
|
np.repeat(by[selected], counts),
|
|
np.repeat(bz[selected], counts),
|
|
)
|
|
cx_r, cy_r, cz_r = (
|
|
np.repeat(cx[selected], counts),
|
|
np.repeat(cy[selected], counts),
|
|
np.repeat(cz[selected], counts),
|
|
)
|
|
denom_r = np.repeat(denom[selected], counts)
|
|
|
|
weight_a = ((by_r - cy_r) * (point_x - cx_r) + (cx_r - bx_r) * (point_y - cy_r)) / denom_r
|
|
weight_b = ((cy_r - ay_r) * (point_x - cx_r) + (ax_r - cx_r) * (point_y - cy_r)) / denom_r
|
|
weight_c = 1.0 - weight_a - weight_b
|
|
hit = (
|
|
(weight_a >= -_BARYCENTRIC_EPS)
|
|
& (weight_b >= -_BARYCENTRIC_EPS)
|
|
& (weight_c >= -_BARYCENTRIC_EPS)
|
|
)
|
|
if not hit.any():
|
|
continue
|
|
z_grid[cell_row[hit], cell_col[hit]] = (
|
|
weight_a[hit] * az_r[hit] + weight_b[hit] * bz_r[hit] + weight_c[hit] * cz_r[hit]
|
|
)
|
|
|
|
return z_grid
|
|
|
|
|
|
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]]:
|
|
"""표현별 npz 모델에서 표고 격자를 환원한 뒤 등고선 리스트를 추출한다."""
|
|
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)
|
|
|
|
if representation == "regular_grid":
|
|
x_coords, y_coords, z_grid, valid_mask = (
|
|
data["x"],
|
|
data["y"],
|
|
data["z"],
|
|
data["valid_mask"],
|
|
)
|
|
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:
|
|
return extract_contours_from_grid(
|
|
x_coords[::step],
|
|
y_coords[::step],
|
|
z_grid[::step, ::step],
|
|
valid_mask[::step, ::step],
|
|
interval,
|
|
scene_center=scene_center,
|
|
)
|
|
return extract_contours_from_grid(
|
|
x_coords, y_coords, z_grid, valid_mask, interval, scene_center=scene_center
|
|
)
|
|
|
|
if representation == "triangular_mesh":
|
|
vertices, faces = data["vertices"], data["faces"]
|
|
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]))
|
|
x_coords, y_coords = _grid_axes(x_min, x_max, y_min, y_max, target_grid_m)
|
|
xx, yy = np.meshgrid(x_coords, y_coords)
|
|
# 저장된 faces를 그대로 쓴다 — griddata는 정점을 다시 삼각분할하느라 긴 변
|
|
# 제거·외곽 클리핑이 반영되지 않은 면까지 되살렸다 (PLAN 2026-08-17).
|
|
z_grid = rasterize_triangle_mesh(vertices, faces, x_coords, y_coords)
|
|
face_mask = _tin_face_coverage_mask(vertices, faces, xx, yy)
|
|
valid_mask = np.isfinite(z_grid) & face_mask
|
|
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
|
|
)
|
|
|
|
if representation == "bspline_surface":
|
|
# 모델 빌더와 같은 곡면을 다시 만든다 — 평활 계수·발산 방어를 공유한다.
|
|
from B04_PreProcess.B04_PreProcess_Engine_ModelContext import (
|
|
evaluate_nurbs_spline,
|
|
fit_nurbs_spline,
|
|
)
|
|
|
|
control_x, control_y, control_z = data["control_x"], data["control_y"], data["control_z"]
|
|
degree = int(data["degree"][0])
|
|
spline, z_range = fit_nurbs_spline(control_x, control_y, control_z, degree)
|
|
x_coords, y_coords = _grid_axes(
|
|
float(control_x[0]),
|
|
float(control_x[-1]),
|
|
float(control_y[0]),
|
|
float(control_y[-1]),
|
|
target_grid_m,
|
|
)
|
|
z_grid = evaluate_nurbs_spline(
|
|
spline, y_coords, x_coords, z_range, Path(model_npz_path).stem
|
|
).astype(np.float32)
|
|
valid_mask = _apply_footprint(
|
|
model_npz_path, x_coords, y_coords, np.ones_like(z_grid, dtype=bool)
|
|
)
|
|
return extract_contours_from_grid(
|
|
x_coords, y_coords, z_grid, valid_mask, interval, scene_center=scene_center
|
|
)
|
|
|
|
if representation == "local_rbf_height_field":
|
|
centers_xy, center_z = data["centers_xy"], 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]))
|
|
x_coords, y_coords = _grid_axes(x_min, x_max, y_min, y_max, target_grid_m)
|
|
xx, yy = np.meshgrid(x_coords, y_coords)
|
|
z_values = interpolator(np.column_stack([xx.ravel(), yy.ravel()])).astype(np.float32)
|
|
z_grid = z_values.reshape(len(y_coords), len(x_coords))
|
|
valid_mask = _apply_footprint(
|
|
model_npz_path, x_coords, y_coords, np.ones_like(z_grid, dtype=bool)
|
|
)
|
|
return extract_contours_from_grid(
|
|
x_coords, y_coords, z_grid, valid_mask, interval, scene_center=scene_center
|
|
)
|
|
|
|
if representation == "meshfree_surfels":
|
|
from scipy.spatial import Delaunay
|
|
|
|
points = data["points"]
|
|
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]))
|
|
x_coords, y_coords = _grid_axes(x_min, x_max, y_min, y_max, target_grid_m)
|
|
# Delaunay 1회로 삼각분할과 hull 판정을 함께 해결한다. 기존에는 griddata가
|
|
# 내부에서 한 번, hull 판정용으로 또 한 번 만들었다 (PLAN 2026-08-17).
|
|
triangulation = Delaunay(np.asarray(points[:, :2], dtype=np.float64))
|
|
z_grid = rasterize_triangle_mesh(points, triangulation.simplices, x_coords, y_coords)
|
|
# hull 밖은 어느 삼각형에도 안 닿아 NaN으로 남는다 = 기존 hull_inside와 같다.
|
|
valid_mask = np.isfinite(z_grid)
|
|
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
|
|
)
|
|
|
|
raise ValueError(f"지원하지 않는 표현 방식입니다: {representation}")
|