diff --git a/B04_PreProcess/B04_PreProcess_Engine_Contour.py b/B04_PreProcess/B04_PreProcess_Engine_Contour.py index 2e793353..3ed70a52 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_Contour.py +++ b/B04_PreProcess/B04_PreProcess_Engine_Contour.py @@ -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