perf(배수유역): 342-3 빠른 길 — 흐름 방향 큰 격자 한 번 · 도로 닿음 판정 한꺼번에

- 하강 방향장 = 해석 격자 둘레 400 m 큰 격자로 한 번 짓고 회차마다 잘라 씀(crop_descent) · 넘으면 그때 다시 지음
- classify_flow 한 칸씩 걷기 → 포인터 더블링(_reach_by_doubling) · 기록용 출발 셀 수는 최외곽 · 나머지 셀 수로
- 거울 시험 test_342_3_watershed_fast(옛 걷기와 칸마다 같음 · 잘라 씀 안쪽 같음) · 실측 bench_342_watershed
- 검증_sub3 실측: 적색 차이 0칸 · 도로 귀속 · 강도 · 관 자리 같음 · 108 → 27 s(붐빈 랩탑 · 340 기준 약 15 s)

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01XWh5wNNycMkbwMpACz1rte
This commit is contained in:
2026-10-16 11:01:28 +09:00
co-authored by Claude Opus 5.5
parent 2879f9bea9
commit 0efc77380e
4 changed files with 349 additions and 65 deletions
@@ -10,6 +10,9 @@
④ **새로 추가된 셀에 적색이 하나도 없으면 종료** (2026-07-31 사용자 지시)
격자 bbox에 닿으면 격자 자체도 셀 정수배로 넓힌다 — 도로 시작점 기준 격자점은 유지된다.
하강 방향장은 처음에 **넉넉한 큰 격자**로 한 번 짓고, 회차마다 그 격자에서 잘라 쓴다
(`DESCENT_SPARE_M` · 342-3). 격자가 커질 때마다 다시 짓던 것(8회 · 30 s)을 한 번으로 줄인다.
"""
from __future__ import annotations
@@ -43,6 +46,11 @@ from config.config_system import (
logger = logging.getLogger(__name__)
# 하강 방향장을 지을 때 해석 격자 둘레에 더 두는 여유(m). 확장이 이 안에서 끝나면 방향장을
# 한 번만 짓는다 — 검증_sub3 실측 확장은 첫 격자에서 남 350 m · 서 200 m 까지 넓어졌다.
# 넘으면 그때 넓힌 격자 둘레로 다시 짓는다(결과는 같고 시간만 듦).
DESCENT_SPARE_M = 400.0
@dataclass
class GridAnalysis:
@@ -78,16 +86,23 @@ def analyze_domain(
"""주어진 격자·해석 영역에 대해 등고선 하강 방향과 도로 도달 색을 한 번 계산한다.
**하강 방향장은 해석 영역과 무관하다** — 등고선 기하만으로 정해진다. 그래서 확장
회차마다 다시 계산하지 않고, 격자가 커졌을 때만 새로 만들어 넘겨받는다(`descent`).
회차마다 다시 계산하지 않고, 이 격자를 덮는 큰 방향장(`descent`)을 받으면 잘라 쓴다.
해석 영역은 마지막에 마스크로만 씌운다.
"""
from B03_FileInput.B03_FileInput_Service_Chain import _log_steps
marks = [("시작", time.perf_counter())]
reused = descent is not None and descent.spec == spec
if not reused:
fitted = crop_descent(descent, spec) if descent is not None else None
if fitted is None:
descent = build_contour_descent(spec, contour_features, None, elevation_floor_m)
marks.append(("하강 방향장" + ("(재사용)" if reused else "(새로 만듦)"), time.perf_counter()))
else:
descent = fitted
marks.append(
(
"하강 방향장" + ("(잘라 씀)" if fitted is not None else "(새로 만듦)"),
time.perf_counter(),
)
)
valid = descent.valid & domain
if not valid.any():
return None
@@ -128,8 +143,12 @@ def expand_by_red_boundary(
# 두면 방향장을 다시 만들지 않고 해석 영역만 넓히며 몇 회차를 돌 수 있다.
spec, domain = _pad_spec(spec, domain, band_cells * 2)
started_cells = int(domain.sum())
spare_cells = max(0, int(round(DESCENT_SPARE_M / spec.cell_m)))
source = build_contour_descent(
_pad_spec(spec, domain, spare_cells)[0], contour_features, None, elevation_floor_m
)
analysis = analyze_domain(
spec, domain, contour_features, route_line, upstream_streams, elevation_floor_m
spec, domain, contour_features, route_line, upstream_streams, elevation_floor_m, source
)
if analysis is None:
return None
@@ -160,6 +179,14 @@ def expand_by_red_boundary(
if not added_mask.any():
closed = True
break
if crop_descent(source, grown_spec) is None:
# 넉넉히 둔 방향장 밖으로 넓어짐 — 넓힌 격자 둘레로 다시 짓는다.
source = build_contour_descent(
_pad_spec(grown_spec, grown_domain, spare_cells)[0],
contour_features,
None,
elevation_floor_m,
)
logger.info(
"배수유역: %d회차 확장 — 최외곽 적색 %d셀 주변 %.0fm, 셀 %d개 추가",
@@ -175,22 +202,20 @@ def expand_by_red_boundary(
route_line,
upstream_streams,
elevation_floor_m,
# 격자가 그대로면 방향장을 재사용한다 — 등고선 기하가 안 바뀌었으므로 결과는 같다.
descent=analysis.descent if grown_spec == current else None,
descent=source,
)
if widened_analysis is None:
break
analysis = widened_analysis
rounds += 1
logger.info(
"[계측] 확장 %d회차 %.1fs · 격자 %d×%d · 셀 %d (+%d) · 방향장 %s",
"[계측] 확장 %d회차 %.1fs · 격자 %d×%d · 셀 %d (+%d)",
attempt + 1,
time.perf_counter() - _round_started,
grown_spec.n_rows,
grown_spec.n_cols,
int(widened.sum()),
int(added_mask.sum()),
"재사용" if grown_spec == current else "새로",
)
added_reaches = widened_analysis.flow.reaches_road.reshape(
@@ -230,6 +255,59 @@ def _dilate_by(mask: np.ndarray, steps: int) -> np.ndarray:
return result
def crop_descent(source: ContourDescent, spec: GridSpec) -> ContourDescent | None:
"""큰 방향장에서 `spec` 격자 몫을 잘라 낸다. 덮지 못하거나 셀이 안 맞으면 None.
잘린 가장자리 밖을 가리키던 셀은 제자리(싱크)로 둔다 — 해석 영역 밖으로 빠지는 물과
같은 판정(미도달)이다. 화살표(방위)는 원래 방향 그대로 둔다.
"""
big = source.spec
if big == spec:
return source
if big.cell_m != spec.cell_m:
return None
row_shift = (big.y_max - spec.y_max) / spec.cell_m
col_shift = (spec.x_min - big.x_min) / spec.cell_m
row0, col0 = int(round(row_shift)), int(round(col_shift))
if (
abs(row_shift - row0) > 1e-6
or abs(col_shift - col0) > 1e-6
or row0 < 0
or col0 < 0
or row0 + spec.n_rows > big.n_rows
or col0 + spec.n_cols > big.n_cols
):
return None
rows = slice(row0, row0 + spec.n_rows)
cols = slice(col0, col0 + spec.n_cols)
def window(flat: np.ndarray) -> np.ndarray:
return flat.reshape(big.n_rows, big.n_cols)[rows, cols]
target = window(source.receiver).astype(np.int64)
target_row = target // big.n_cols - row0
target_col = target % big.n_cols - col0
inside = (
(target_row >= 0)
& (target_row < spec.n_rows)
& (target_col >= 0)
& (target_col < spec.n_cols)
)
local = np.arange(spec.size, dtype=np.int64).reshape(spec.n_rows, spec.n_cols)
receiver = np.where(inside, target_row * spec.n_cols + target_col, local)
return ContourDescent(
spec=spec,
band_elevation=source.band_elevation[rows, cols].copy(),
valid=source.valid[rows, cols].copy(),
receiver=receiver.reshape(-1).astype(np.int32),
step_length=np.where(inside, window(source.step_length), 0.0)
.reshape(-1)
.astype(np.float32),
azimuth=window(source.azimuth).reshape(-1).copy(),
levels=source.levels,
)
def _pad_spec(spec: GridSpec, domain: np.ndarray, cells: int) -> tuple[GridSpec, np.ndarray]:
"""격자에 사방 여유를 두고 해석 영역 마스크를 그 안으로 옮겨 담는다."""
if cells <= 0:
@@ -278,8 +278,8 @@ class FlowClassification:
direction: np.ndarray # (R*C,) int16 — 32방위 코드(32=제자리, 33=무효)
reaches_road: np.ndarray # (R*C,) bool — 물길을 따라가면 도로에 닿는가
analyzed: np.ndarray # (R*C,) bool — 실제로 판정한 셀
outer_seeds: int # 최외곽에서 출발해 판정한 셀 수
interior_seeds: int # 최외곽 추적에 안 걸려 따로 출발시킨 내부 셀 수
outer_seeds: int # 최외곽 셀 수(기록용 · 세류 셀 뺌)
interior_seeds: int # 그 밖 판정 셀 수(기록용 · 세류 셀 뺌)
burned: np.ndarray | None = None # (R*C,) bool — 세류망을 따라 흐름을 새긴 셀
@@ -321,6 +321,8 @@ def classify_flow(
다시 분석하지 않는다.
④ 최외곽 추적에 안 걸린 내부 셀을 그다음에 따로 출발시킨다.
판정 값은 이 차례로 걷는 것과 칸마다 같다 — 계산은 한꺼번에 한다(`_reach_by_doubling`).
`azimuth`를 주면 그 32방위 코드를 그대로 화살표로 쓴다(등고선 하강 방향). 주지 않으면
지표면 기울기에서 뽑는다. **화살표는 실제 수신 셀과 같은 방향이어야 한다** — 어긋나면
화살표로 사슬을 따라가는 눈 검증이 성립하지 않는다.
@@ -336,21 +338,13 @@ def classify_flow(
direction = descent_azimuth(spec, terrain.elevation, terrain.valid, receiver, burned)
else:
direction = _azimuth_with_burn(spec, azimuth, receiver, burned)
reaches = np.zeros(spec.size, dtype=bool)
# 0=미방문, 1=경로에 올라 있음, 2=판정 완료
state = np.zeros(spec.size, dtype=np.int8)
# ⓪ 세류선과 겹치는 셀을 먼저 적색으로 못박는다.
reaches[stream_cells] = True
state[stream_cells] = 2
reaches = _reach_by_doubling(receiver, valid, absorbing, stream_cells)
analyzed = valid | stream_cells
stream_done = int(stream_cells.sum())
# 기록용 숫자 — 한꺼번에 판정하므로 「출발」 구분이 없다. 최외곽 셀 수 · 나머지 셀 수만 남긴다.
outer_done = int((outermost_cells(terrain.valid).reshape(-1) & ~stream_cells).sum())
interior_done = int((valid & ~stream_cells).sum()) - outer_done
outer = np.flatnonzero(outermost_cells(terrain.valid))
remaining = np.flatnonzero(valid)
outer_done = _walk_from(outer, receiver, valid, absorbing, state, reaches)
interior_done = _walk_from(remaining, receiver, valid, absorbing, state, reaches)
analyzed = state == 2
logger.info(
"배수유역: 흐름 판정 %d셀 (세류 선확정 %d / 최외곽 출발 %d / 내부 보충 %d) — "
"도로 도달 %d, 미도달 %d",
@@ -391,52 +385,32 @@ def _azimuth_with_burn(
return direction
def _walk_from(
starts: np.ndarray,
def _reach_by_doubling(
receiver: np.ndarray,
valid: np.ndarray,
absorbing: np.ndarray,
state: np.ndarray,
reaches: np.ndarray,
) -> int:
"""출발 셀에서 물길 사슬을 따라가며 판정하고, 사슬 전체에 같은 색을 적는다.
stream_cells: np.ndarray,
) -> np.ndarray:
"""물길 사슬을 끝까지 따라간 종착 셀로 셀마다 도로 도달을 한꺼번에 판정한다.
`absorbing`은 도로 셀과 확정된 세류망 셀 — 여기 닿으면 사슬 전체가 도로 도달이다.
새로 판정한 셀 수를 돌려준다.
사슬을 한 칸씩 걷던 판정(최외곽 출발 → 내부 보충 · 이미 판정된 셀의 색 물려받기)과
**칸마다 같은 값**이다 — 순환이 없으면 어느 셀에서 출발하든 같은 종착에 닿기 때문이다.
종착 = 도로 · 세류 셀(적색) · 싱크 · 해석 영역 밖으로 나가는 셀(미도달). 세류 셀은 표고가
없어도 따라간다. 포인터 더블링(`trace_flow` 와 같은 방식)이라 파이썬 반복이 없다.
"""
resolved = 0
path: list[int] = []
for start in starts.tolist():
if state[start] == 2:
continue
path.clear()
node = start
while True:
if state[node] == 2:
verdict = bool(reaches[node]) # 이미 색이 정해진 셀 — 그 색을 물려받는다
break
if state[node] == 1: # 방어: 채움·평탄해소 후에는 순환이 없어야 한다
verdict = False
break
state[node] = 1
path.append(node)
if absorbing[node]:
verdict = True # 도로 또는 세류망에 합류 — 여기서 하류로 빠진다
break
following = int(receiver[node])
if following == node:
verdict = False # 싱크에 갇힘
break
if not valid[following] and state[following] != 2:
# 해석 영역 밖으로 빠짐. 단, 이미 색이 정해진 셀(세류 선확정 등)이면 따라간다.
verdict = False
break
node = following
for visited in path:
reaches[visited] = verdict
state[visited] = 2
resolved += len(path)
return resolved
index = np.arange(receiver.size, dtype=np.int64)
following = receiver.astype(np.int64)
terminal = absorbing | (following == index) | ~(valid[following] | stream_cells[following])
jump = np.where(terminal, index, following)
for _ in range(_MAX_DOUBLING_ROUNDS):
next_jump = jump[jump]
if np.array_equal(next_jump, jump):
break
jump = next_jump
# 순환에 갇힌 사슬은 종착에 못 닿는다 — 걷던 판정처럼 미도달.
reaches = absorbing[jump] & terminal[jump] & (valid | stream_cells)
reaches[stream_cells] = True
return reaches
# ── 격자 확장 (별도 단계) ───────────────────────────────────────────────────
+94
View File
@@ -0,0 +1,94 @@
"""342-3 · 347-1 배수유역 빠른 길 실측 — 검증_sub3 도엽 · 노선으로 preview_stages 를 돌려
시간과 결과(적색 · 흐름 방향 · 도로 귀속 · 강도)를 tmp/342/<이름>.npz 에 남기고 기준과 대조한다.
./venv/Scripts/python.exe resources/tester/bench_342_watershed.py base # 기준 저장
./venv/Scripts/python.exe resources/tester/bench_342_watershed.py new base # 새로 돌려 대조
pytest 대상 아님(DB · 프로젝트 파일이 있어야 돎). 저장 안 함 — 프로젝트 파일은 읽기만
(흐름 방향 저장본 drainage/descent_*.npz 는 347-1 부터 씀).
"""
import asyncio
import json
import logging
import sys
import time
from pathlib import Path
from uuid import UUID
import numpy as np
ROOT = Path(__file__).resolve().parents[2]
sys.path.insert(0, str(ROOT))
PID = UUID("1f93a515-df0a-4975-a8ba-daaa7672b123") # 검증_sub3
OUT_DIR = ROOT / "tmp" / "342"
async def _prepare():
from config.config_db import close_db_pool, init_db_pool
await init_db_pool()
try:
from B04_PreProcess.B04_PreProcess_Router_Watershed import _prepare as prepare
started = time.perf_counter()
prepared = await prepare(PID)
return prepared, time.perf_counter() - started
finally:
await close_db_pool()
def main() -> None:
name = sys.argv[1]
base = sys.argv[2] if len(sys.argv) > 2 else None
OUT_DIR.mkdir(parents=True, exist_ok=True)
logging.basicConfig(
level=logging.INFO,
format="%(asctime)s %(name)s:%(message)s",
handlers=[logging.FileHandler(OUT_DIR / f"{name}.log", "w", encoding="utf-8")],
)
prepared, prepare_s = asyncio.run(_prepare())
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Analyze import preview_stages
started = time.perf_counter()
pv = preview_stages(prepared["vertices"], prepared["contours"], prepared["streams"])
preview_s = time.perf_counter() - started
spec = pv.spec
red = pv.flow.reaches_road & pv.flow.analyzed
result = {
"spec": np.array([spec.x_min, spec.y_max, spec.cell_m, spec.n_rows, spec.n_cols]),
"domain": pv.domain,
"red": red,
"direction": pv.flow.direction,
"road_slot": pv.routing.road_slot,
"strength": np.array(pv.strength_profile),
"pipes": np.array([p.chainage_m for p in pv.pipes]),
}
np.savez_compressed(OUT_DIR / f"{name}.npz", **result)
report = {
"prepare_s": round(prepare_s, 2),
"preview_s": round(preview_s, 2),
"total_s": round(prepare_s + preview_s, 2),
"rounds": pv.expand_rounds,
"grid": [spec.n_rows, spec.n_cols],
"red_cells": int(red.sum()),
}
if base:
ref = np.load(OUT_DIR / f"{base}.npz")
same_spec = bool(np.array_equal(ref["spec"], result["spec"]))
report["same_spec"] = same_spec
if same_spec:
for key in ("domain", "red", "direction", "road_slot"):
report[f"diff_{key}"] = int((ref[key] != result[key]).sum())
else:
report["base_spec"] = ref["spec"].tolist()
report["new_spec"] = result["spec"].tolist()
report["red_cells_base"] = int(ref["red"].sum())
report["same_strength"] = bool(np.array_equal(ref["strength"], result["strength"]))
report["same_pipes"] = bool(np.array_equal(ref["pipes"], result["pipes"]))
print(json.dumps(report, ensure_ascii=False))
(OUT_DIR / f"{name}.json").write_text(json.dumps(report, ensure_ascii=False, indent=1), "utf-8")
if __name__ == "__main__":
main()
@@ -0,0 +1,138 @@
"""342-3 배수유역 빠른 길 — 도로 닿음 판정 한꺼번에 · 큰 방향장 잘라 씀.
- 판정: 옛 한 칸씩 걷기(최외곽 출발 → 내부 보충 · 색 물려받기)와 칸마다 같음(적색 차이 0).
- 잘라 씀: 같은 격자면 그대로 · 셀 어긋남 · 범위 밖은 None · 안쪽은 제 격자로 지은 것과 같음.
실측(검증_sub3 · 적색 차이 0 · 60 → 약 15 s)은 bench_342_watershed.py.
"""
import math
import numpy as np
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Descent import build_contour_descent
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Expand import _pad_spec, crop_descent
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Flow import (
_reach_by_doubling,
outermost_cells,
)
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import GridSpec
def _walk_reference(receiver, valid2d, absorbing, stream):
"""옛 classify_flow 판정(한 칸씩 걷기) 그대로."""
valid = valid2d.reshape(-1)
reaches = np.zeros(receiver.size, bool)
state = np.zeros(receiver.size, np.int8)
reaches[stream] = True
state[stream] = 2
starts = list(np.flatnonzero(outermost_cells(valid2d))) + list(np.flatnonzero(valid))
for start in starts:
if state[start] == 2:
continue
path, node = [], start
while True:
if state[node] == 2:
verdict = bool(reaches[node])
break
if state[node] == 1:
verdict = False
break
state[node] = 1
path.append(node)
if absorbing[node]:
verdict = True
break
following = int(receiver[node])
if following == node:
verdict = False
break
if not valid[following] and state[following] != 2:
verdict = False
break
node = following
for visited in path:
reaches[visited] = verdict
state[visited] = 2
return reaches
def _random_case(seed, rows=60, cols=70):
rng = np.random.default_rng(seed)
size = rows * cols
potential = rng.random(size)
index = np.arange(size)
r, c = index // cols, index % cols
receiver = index.copy()
best = potential.copy()
for dr in (-1, 0, 1):
for dc in (-1, 0, 1):
nr, nc = r + dr, c + dc
ok = (nr >= 0) & (nr < rows) & (nc >= 0) & (nc < cols)
n = np.where(ok, nr * cols + nc, index)
better = ok & (potential[n] < best)
best = np.where(better, potential[n], best)
receiver = np.where(better, n, receiver)
valid2d = rng.random((rows, cols)) > 0.15
road = np.zeros(size, bool)
road[(r == rows // 2) & (c > 5) & (c < cols - 5)] = True
stream = (rng.random(size) > 0.97) & ~road # 표고 없는 세류 셀도 섞임
# 세류 셀은 물길 다음 셀로 새김 — 순환 하나를 일부러 만든다(걷던 판정 = 미도달).
receiver = receiver.astype(np.int32)
receiver[stream] = np.roll(index, 1)[stream]
receiver[5 * cols + 5], receiver[5 * cols + 6] = 5 * cols + 6, 5 * cols + 5
absorbing = (road & valid2d.reshape(-1)) | stream
return receiver, valid2d, absorbing, stream
def test_reach_matches_walk_cell_by_cell():
for seed in range(8):
receiver, valid2d, absorbing, stream = _random_case(seed)
expected = _walk_reference(receiver, valid2d, absorbing, stream)
got = _reach_by_doubling(receiver, valid2d.reshape(-1), absorbing, stream)
assert int((got != expected).sum()) == 0, seed
assert expected.sum() > 0 and (~expected).sum() > 0
def _circle_contours(cx, cy, radii_levels):
features = []
for radius, level in radii_levels:
coords = [
[cx + radius * math.cos(t), cy + radius * math.sin(t)]
for t in np.linspace(0, 2 * math.pi, 181)
]
features.append(
{
"geometry": {"type": "LineString", "coordinates": coords},
"properties": {"elevation": level},
}
)
return features
def test_crop_descent_matches_own_grid_inside():
contours = _circle_contours(100.0, -100.0, [(r, 200.0 - r) for r in range(10, 140, 10)])
spec = GridSpec(x_min=40.0, y_max=-40.0, cell_m=1.0, n_rows=120, n_cols=120)
big_spec, _ = _pad_spec(spec, np.ones((120, 120), bool), 30)
big = build_contour_descent(big_spec, contours)
own = build_contour_descent(spec, contours)
cut = crop_descent(big, spec)
assert cut is not None and cut.spec == spec
inner = np.zeros((120, 120), bool)
inner[20:100, 20:100] = True
flat = inner.reshape(-1)
assert np.array_equal(cut.valid[inner], own.valid[inner])
assert np.array_equal(cut.receiver[flat], own.receiver[flat])
assert np.array_equal(cut.azimuth[flat], own.azimuth[flat])
# 잘린 몫 밖을 가리키는 셀은 없다(제자리로 둠).
assert cut.receiver.min() >= 0 and cut.receiver.max() < spec.size
assert crop_descent(big, big_spec) is big
def test_crop_descent_refuses_misaligned_or_outside():
contours = _circle_contours(100.0, -100.0, [(r, 200.0 - r) for r in range(10, 140, 10)])
spec = GridSpec(x_min=40.0, y_max=-40.0, cell_m=1.0, n_rows=50, n_cols=50)
big = build_contour_descent(spec, contours)
shifted = GridSpec(x_min=40.5, y_max=-40.0, cell_m=1.0, n_rows=20, n_cols=20)
outside = GridSpec(x_min=30.0, y_max=-40.0, cell_m=1.0, n_rows=20, n_cols=20)
assert crop_descent(big, shifted) is None
assert crop_descent(big, outside) is None