eomsangdon/Aislo
1
0
Fork 0
Code Issues Pull Requests 2 Actions Packages Projects Releases Wiki Activity

perf(B04): 등고선 추출 삼각망 래스터화 — 15유닛 27분39초에서 77.6초로

Browse Source 
griddata가 질의점마다 삼각형 98만개 망에서 위치탐색(find_simplex)을 돌던
것이 등고선 1회 2~5분의 원인이었다. 삼각형을 격자에 직접 래스터화해
위치탐색을 제거한다. 분석 범위 15유닛과 스무딩 온오프는 그대로 유지.

- rasterize_triangle_mesh 신설: 삼각형별 격자 bbox만 순회하며 무게중심
  선형보간. 격자 인덱스는 np.searchsorted로 축을 직접 탐색한다 — _grid_axes가
  float32 축을 만드는데 UTM 좌표에서 해상도가 0.015625m라 등간격 환산 시
  최대 0.81셀 어긋나 실데이터 8,819셀을 놓쳤다.
- triangular_mesh: griddata 제거, npz에 저장된 faces 재사용. 저장 faces는
  긴 변 제거·외곽 클리핑이 반영된 실제 TIN이라 재삼각분할보다 충실하다.
- meshfree_surfels: Delaunay 1회로 삼각분할과 hull 판정을 함께 해결.
  기존의 griddata 내부 + 별도 Delaunay 2회를 1회로 줄였다.

검증: meshfree·dtm 산출물 완전 일치(세그먼트·정점 수), tin은 삼각망 차이만큼
정점 -0.03%(레벨 집합·좌표 범위 동일, 표고차 최대 0.4083m). 래스터화에 같은
재Delaunay를 넣으면 griddata와 셀 단위 차이 0 — 보간 정확성 확인.
tmp/tests 66건 통과.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This commit is contained in:
eomsangdon
2026-08-17 02:03:03 +09:00
co-authored by Claude Opus 5
parent 84eb027348
commit 5d0346aa5c
1 changed files with 135 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files
B04_PreProcess/B04_PreProcess_Engine_Contour.py
+135 -12
View File
@@ -219,6 +219,133 @@ def _grid_axes(x_min: float, x_max: float, y_min: float, y_max: float, target_gr
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,
@@ -255,14 +382,14 @@ def extract_contours(
)
if representation == "triangular_mesh":
from scipy.interpolate import griddata
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)
z_grid = griddata(vertices[:, :2], vertices[:, 2], (xx, yy), method="linear")
# 저장된 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)
@@ -320,22 +447,18 @@ def extract_contours(
)
if representation == "meshfree_surfels":
from scipy.interpolate import griddata
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)
xx, yy = np.meshgrid(x_coords, y_coords)
z_grid = griddata(points[:, :2], points[:, 2], (xx, yy), method="linear")
# 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)
try:
tri = Delaunay(points[:, :2])
hull_inside = tri.find_simplex(np.column_stack([xx.ravel(), yy.ravel()])) >= 0
valid_mask = valid_mask & hull_inside.reshape(xx.shape)
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