diff --git a/B04_PreProcess/B04_PreProcess_Engine_Watershed_Descent.py b/B04_PreProcess/B04_PreProcess_Engine_Watershed_Descent.py index 14a00173..b15c16ca 100644 --- a/B04_PreProcess/B04_PreProcess_Engine_Watershed_Descent.py +++ b/B04_PreProcess/B04_PreProcess_Engine_Watershed_Descent.py @@ -121,6 +121,90 @@ def rasterize_contours( return burned, levels +# 밴드별 하강거리 EDT 를 **쓸 자리 둘레**에서만 돌리기 위한 값들. +# +# 왜(2026-09-07 실측) — 확장 회차마다 방향장을 다시 만드는데, 그 안의 밴드별 EDT 가 +# 격자 전체를 단계마다 훑어 확장 95.8s 중 50.4s 를 썼다. 정작 쓰는 값은 그 단계의 +# 밴드 셀뿐이다. 다만 `analyze_domain` 이 `domain` 을 안 넘겨 밴드가 **격자 전역에 +# 흩어지므로** 바운딩박스 하나로는 안 좁아진다. 그래서 격자를 타일로 나눠 **셀이 있는 +# 타일만** 그 둘레 `margin` 까지 잘라 EDT 를 돌린다. +# +# 창 밖에 더 가까운 등고선이 있을 수 있으면(창 안 최대 거리가 여유에 닿거나 창에 낮은 +# 라인이 없으면) 그 단계를 전체 격자로 다시 돌린다 — **근사가 아니라 같은 값을 싸게 +# 구하는 것**. 값이 같은지는 `np.array_equal` 로 확인했다(거리·목표행·목표열 전부). +# +# 고른 근거(745×1035 격자·166단 중 EDT 도는 91단, 전체격자 방식 10.8s): +# 타일 128 여유 32 → 3.2s(폴백 9) 타일 256 여유 24 → 3.3s(폴백 27) +# **타일 256 여유 32 → 2.9s(폴백 9)** 타일 256 여유 48 → 3.2s(폴백 4) +# 타일 512 여유 32 → 4.4s(폴백 9) +# 여유는 타일 크기와 무관하게 **거리 실측**으로 정했다 — 한 단 아래 등고선은 주곡선 +# 5m·셀 1m 에서 대개 수십 칸 안이다. 지형이 완만해 폴백이 잦아지면(로그에 횟수가 +# 찍힌다) 여유를 늘릴 것. +DESCENT_WINDOW_MARGIN_CELLS = 32 +DESCENT_TILE_CELLS = 256 + + +def _distance_to_lower( + members: np.ndarray, + lower: np.ndarray, + margin: int | None = None, + tile: int | None = None, +) -> tuple[np.ndarray, np.ndarray, np.ndarray, bool]: + """`members` 셀에서 한 단 낮은 등고 라인까지의 거리와 목표 셀(**전체 격자 좌표**). + + 네 번째 값은 전체 격자로 되돌아갔는지(폴백) 여부다. 반환 순서는 `members` 의 + 행우선 순서 — 호출부가 `distance[members] = ...` 로 그대로 넣는다. + """ + margin = DESCENT_WINDOW_MARGIN_CELLS if margin is None else margin + tile = DESCENT_TILE_CELLS if tile is None else tile + rows, cols = members.shape + out_distance = np.zeros((rows, cols), dtype=np.float64) + out_row = np.zeros((rows, cols), dtype=np.int32) + out_col = np.zeros((rows, cols), dtype=np.int32) + + for row_start in range(0, rows, tile): + row_stop = min(rows, row_start + tile) + for col_start in range(0, cols, tile): + col_stop = min(cols, col_start + tile) + tile_members = members[row_start:row_stop, col_start:col_stop] + if not tile_members.any(): + continue + win_row0 = max(0, row_start - margin) + win_row1 = min(rows, row_stop + margin) + win_col0 = max(0, col_start - margin) + win_col1 = min(cols, col_stop + margin) + window_lower = lower[win_row0:win_row1, win_col0:win_col1] + if not window_lower.any(): + return _full_distance_to_lower(members, lower) + step, (step_row, step_col) = distance_transform_edt(~window_lower, return_indices=True) + in_window = np.zeros(window_lower.shape, dtype=bool) + in_window[ + row_start - win_row0 : row_stop - win_row0, + col_start - win_col0 : col_stop - win_col0, + ] = tile_members + hit_rows, hit_cols = np.nonzero(in_window) + values = step[hit_rows, hit_cols] + # 타일 둘레로 `margin` 을 뒀으므로, 거리가 그보다 짧으면 창 밖에 더 가까운 + # 것은 있을 수 없다. 닿으면 이 단계는 통째로 전체 격자로 다시 잰다. + if values.size and float(values.max()) >= margin: + return _full_distance_to_lower(members, lower) + global_rows = hit_rows + win_row0 + global_cols = hit_cols + win_col0 + out_distance[global_rows, global_cols] = values + out_row[global_rows, global_cols] = step_row[hit_rows, hit_cols] + win_row0 + out_col[global_rows, global_cols] = step_col[hit_rows, hit_cols] + win_col0 + + return out_distance[members], out_row[members], out_col[members], False + + +def _full_distance_to_lower( + members: np.ndarray, lower: np.ndarray +) -> tuple[np.ndarray, np.ndarray, np.ndarray, bool]: + """격자 전체 EDT — 창으로 못 믿을 때만 쓰는 폴백.""" + step, (step_row, step_col) = distance_transform_edt(~lower, return_indices=True) + return step[members], step_row[members], step_col[members], True + + def build_contour_descent( spec: GridSpec, contour_features: list[dict[str, Any]], @@ -160,6 +244,7 @@ def build_contour_descent( target_row = np.zeros((rows, cols), dtype=np.int32) target_col = np.zeros((rows, cols), dtype=np.int32) band_rank = np.full((rows, cols), -1, dtype=np.int32) + window_fallbacks = 0 for rank, elevation in enumerate(levels[:-1]): members = inside & (band_elevation == elevation) if not members.any(): @@ -167,15 +252,21 @@ def build_contour_descent( lower = on_contour & (burned < elevation) if not lower.any(): continue - step_distance, (step_row, step_col) = distance_transform_edt(~lower, return_indices=True) - distance[members] = step_distance[members].astype(np.float32) - target_row[members] = step_row[members] - target_col[members] = step_col[members] + step_distance, step_row, step_col, fell_back = _distance_to_lower(members, lower) + window_fallbacks += int(fell_back) + distance[members] = step_distance.astype(np.float32) + target_row[members] = step_row + target_col[members] = step_col band_rank[members] = len(levels) - 1 - rank # 높을수록 큰 값 # 최하단 밴드는 더 내려갈 등고 라인이 없다. 무효로 빼지 않고 **정지 셀**로 남긴다 — # 그래야 그 자리가 도로·세류선이면 적색으로 잡히고, 아니면 물이 고이는 지점으로 보인다. - marks.append((f"밴드별 하강거리(EDT {len(levels) - 1}회)", time.perf_counter())) + marks.append( + ( + f"밴드별 하강거리(EDT {len(levels) - 1}회 · 전체격자 폴백 {window_fallbacks}회)", + time.perf_counter(), + ) + ) lowest = inside & (band_elevation == levels[-1]) & (band_rank < 0) if lowest.any(): distance[lowest] = 0.0