perf(공용 지형): 라플라스 홀짝 네 덩이 · 거리 비례 copyto · 폐합 채움 수위 상자 창 — 저장본과 비트 동일 (PLAN 617-2 · 요구 504 · 506)
- 새 common_util_terrain_laplace.py — B04 build_laplace 짝(같은 덧셈 차례 · float32) · build_contour_surface 가 부름 - 폐합 채움 = 수위 칸 상자 + 1 칸 덧댐 안 채움 · 묶음 find_objects 창 · 차례 그대로 · 쓰지 않는 _window 지움 - TERRAIN_ALGO 1 그대로 · B 파일 안 고침 - 짝 시험 test_617_2_laplace_mirror.py Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Claude Opus 5.5
@@ -0,0 +1,111 @@
|
||||
"""등고 격자 → 라플라스 면 — 공용 지형 창고가 쓰는 빠른 판(PLAN 617-2 · 요구 504 · 506).
|
||||
|
||||
B04 짝: `B04_PreProcess_Engine_SheetMethods.build_laplace`(거리 비례 + `relax_laplace` 400 회)와
|
||||
비트 동일 · 짝 시험 `test_617_2_laplace_mirror.py` · B 파일은 안 고침.
|
||||
빠른 까닭 둘:
|
||||
① 라플라스 — 격자를 홀짝 네 덩이(`S[p::2, q::2]` · 연속 메모리)로 나눠 그 색 칸만 셈 ·
|
||||
덧셈 차례(위 + 아래 + 왼 + 오른) × 0.25 는 원본과 같음.
|
||||
② 거리 비례 — `np.where` 새 배열 대신 `np.copyto(where=)`.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
def distance_surface(burned: np.ndarray, cell_m: float) -> np.ndarray:
|
||||
"""B04 `build_distance` 와 같은 값 — 가장 가까운 서로 다른 두 수위 사이 거리 비례."""
|
||||
from scipy.ndimage import distance_transform_edt
|
||||
|
||||
levels = np.unique(burned[np.isfinite(burned)])
|
||||
if len(levels) < 2:
|
||||
return np.full(burned.shape, levels[0] if len(levels) else np.nan, dtype=np.float32)
|
||||
|
||||
infinity = np.float32(np.inf)
|
||||
first_d = np.full(burned.shape, infinity, dtype=np.float32)
|
||||
first_z = np.zeros(burned.shape, dtype=np.float32)
|
||||
second_d = np.full(burned.shape, infinity, dtype=np.float32)
|
||||
second_z = np.zeros(burned.shape, dtype=np.float32)
|
||||
beats_second = np.empty(burned.shape, dtype=bool)
|
||||
for level in levels:
|
||||
distance = distance_transform_edt(burned != level, sampling=cell_m).astype(np.float32)
|
||||
beats_first = distance < first_d
|
||||
np.copyto(second_d, first_d, where=beats_first)
|
||||
np.copyto(second_z, first_z, where=beats_first)
|
||||
np.copyto(first_d, distance, where=beats_first)
|
||||
np.copyto(first_z, np.float32(level), where=beats_first)
|
||||
np.less(distance, second_d, out=beats_second)
|
||||
beats_second &= ~beats_first
|
||||
np.copyto(second_d, distance, where=beats_second)
|
||||
np.copyto(second_z, np.float32(level), where=beats_second)
|
||||
|
||||
total = first_d + second_d
|
||||
usable = np.isfinite(second_d) & (total > 1e-9)
|
||||
surface = first_z.astype(np.float32)
|
||||
surface[usable] = (
|
||||
(
|
||||
first_z[usable].astype(np.float64) * second_d[usable].astype(np.float64)
|
||||
+ second_z[usable].astype(np.float64) * first_d[usable].astype(np.float64)
|
||||
)
|
||||
/ total[usable].astype(np.float64)
|
||||
).astype(np.float32)
|
||||
return surface
|
||||
|
||||
|
||||
def relax_checker(surface: np.ndarray, fixed: np.ndarray, iterations: int) -> None:
|
||||
"""B04 `relax_laplace` 와 같은 값(in-place · 빨강 → 검정) — 홀짝 네 덩이로 셈.
|
||||
|
||||
덩이 (p, q) 마다 가장자리 한 칸 덧댄 버퍼.
|
||||
(p, q) 칸의 위 · 아래 = (1-p, q) 버퍼의 연이은 두 행 ·
|
||||
왼 · 오른 = (p, 1-q) 버퍼의 연이은 두 열 · 격자 끝은 자기 값(원본 edge 덧대기) → 셈 직전
|
||||
그 덧댄 줄만 자기 첫 · 끝 줄로 채움.
|
||||
"""
|
||||
if iterations <= 0:
|
||||
return
|
||||
free = ~fixed & np.isfinite(surface)
|
||||
if not free.any():
|
||||
return
|
||||
quads = {(p, q): (p + q) & 1 for p in (0, 1) for q in (0, 1)} # 0 = 빨강
|
||||
buffers: dict[tuple[int, int], np.ndarray] = {}
|
||||
for p, q in quads:
|
||||
block = surface[p::2, q::2]
|
||||
buffer = np.zeros((block.shape[0] + 2, block.shape[1] + 2), dtype=np.float32)
|
||||
buffer[1:-1, 1:-1] = block
|
||||
buffers[p, q] = buffer
|
||||
frees = {key: np.ascontiguousarray(free[key[0] :: 2, key[1] :: 2]) for key in quads}
|
||||
temps = {key: np.empty(buffers[key][1:-1, 1:-1].shape, dtype=np.float32) for key in quads}
|
||||
quarter = np.float32(0.25)
|
||||
for _ in range(iterations):
|
||||
for colour in (0, 1):
|
||||
for (p, q), quad_colour in quads.items():
|
||||
if quad_colour != colour:
|
||||
continue
|
||||
own = buffers[p, q][1:-1, 1:-1]
|
||||
rows, cols = own.shape
|
||||
if own.size == 0:
|
||||
continue
|
||||
vertical = buffers[1 - p, q]
|
||||
horizontal = buffers[p, 1 - q]
|
||||
if p == 0:
|
||||
vertical[0, 1:-1] = own[0]
|
||||
if p + rows == vertical.shape[0] - 1:
|
||||
vertical[-1, 1:-1] = own[-1]
|
||||
if q == 0:
|
||||
horizontal[1:-1, 0] = own[:, 0]
|
||||
if q + cols == horizontal.shape[1] - 1:
|
||||
horizontal[1:-1, -1] = own[:, -1]
|
||||
t = temps[p, q]
|
||||
np.add(vertical[p : p + rows, 1:-1], vertical[p + 1 : p + 1 + rows, 1:-1], out=t)
|
||||
t += horizontal[1:-1, q : q + cols]
|
||||
t += horizontal[1:-1, q + 1 : q + 1 + cols]
|
||||
t *= quarter
|
||||
np.copyto(own, t, where=frees[p, q])
|
||||
for p, q in quads:
|
||||
surface[p::2, q::2] = buffers[p, q][1:-1, 1:-1]
|
||||
|
||||
|
||||
def laplace_surface(burned: np.ndarray, cell_m: float, iterations: int = 400) -> np.ndarray:
|
||||
"""B04 `build_laplace` 와 같은 값 — 거리 비례 초기값 + 라플라스 다듬기."""
|
||||
surface = distance_surface(burned, cell_m)
|
||||
relax_checker(surface, np.isfinite(burned), iterations)
|
||||
return surface
|
||||
@@ -7,7 +7,9 @@ B04 `B04_PreProcess_Engine_SheetSurface.py` 라플라스 길을 본뜸 · B 파
|
||||
최대 9 m) → 선분 끝을 정렬한 뒤 큰 축 칸 중심마다 한 칸 + 꼭짓점 칸(방향 무관 · 굵기 1 칸).
|
||||
② 폐합 안쪽 채움 — 봉우리 · 웅덩이 방향 · 기울기를 같은 높이 선 전체가 아니라 그 고리 둘레로만
|
||||
→ 범위가 바뀌어도 같은 고리는 같은 값(조각 = 통째).
|
||||
보간 자체(거리 비례 + 라플라스 400 회)는 B04 `build_laplace` 를 그대로 부름.
|
||||
보간 자체(거리 비례 + 라플라스 400 회)는 `common_util_terrain_laplace`
|
||||
(B04 `build_laplace` 짝 · 비트 동일).
|
||||
폐합 채움은 수위 칸 상자 · 묶음 창 안에서만 셈(통격자 셈과 같은 값 · PLAN 617-2).
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
@@ -17,8 +19,8 @@ from typing import Any
|
||||
|
||||
import numpy as np
|
||||
|
||||
from B04_PreProcess.B04_PreProcess_Engine_SheetMethods import build_laplace
|
||||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import GridSpec
|
||||
from common_util.common_util_terrain_laplace import laplace_surface
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
@@ -90,15 +92,6 @@ def burn_segments(spec: GridSpec, a: np.ndarray, b: np.ndarray, z: np.ndarray) -
|
||||
return burned.astype(np.float32)
|
||||
|
||||
|
||||
def _window(component: np.ndarray, pad: int) -> tuple[slice, slice]:
|
||||
rows = np.flatnonzero(component.any(axis=1))
|
||||
cols = np.flatnonzero(component.any(axis=0))
|
||||
return (
|
||||
slice(max(rows[0] - pad, 0), rows[-1] + pad + 1),
|
||||
slice(max(cols[0] - pad, 0), cols[-1] + pad + 1),
|
||||
)
|
||||
|
||||
|
||||
def resolve_enclosed_interiors(
|
||||
burned: np.ndarray, present: list[float], surface: np.ndarray, cell_m: float, interval_m: float
|
||||
) -> int:
|
||||
@@ -106,27 +99,54 @@ def resolve_enclosed_interiors(
|
||||
|
||||
B04 `_resolve_enclosed_interiors` 와 같은 규칙 · 다른 점 = 방향(바깥 평균)과 기울기 띠를
|
||||
그 칸 묶음을 감싼 고리 칸 바깥으로만 잼(같은 높이 다른 선은 안 섞음). 채운 곳 수를 돌려줌.
|
||||
채움 · 묶음은 그 수위 칸 상자 + 1 칸 빈 덧댐 안에서 · 묶음은 `find_objects` 창 —
|
||||
통격자와 같은 값 · 차례(수위 → 묶음 번호) 그대로.
|
||||
"""
|
||||
from scipy.ndimage import binary_dilation, binary_fill_holes, distance_transform_edt, label
|
||||
from scipy.ndimage import (
|
||||
binary_dilation,
|
||||
binary_fill_holes,
|
||||
distance_transform_edt,
|
||||
find_objects,
|
||||
label,
|
||||
)
|
||||
|
||||
constrained = np.isfinite(burned)
|
||||
limit = max(interval_m - 0.5, 0.5)
|
||||
pad = _BAND_CELLS + 2
|
||||
eight = np.ones((3, 3), dtype=bool)
|
||||
resolved = 0
|
||||
rows_all, cols_all = np.nonzero(constrained)
|
||||
values = burned[rows_all, cols_all]
|
||||
order = np.argsort(values, kind="stable")
|
||||
rows_all, cols_all, values = rows_all[order], cols_all[order], values[order]
|
||||
for level in present:
|
||||
mask = burned == level
|
||||
interior = binary_fill_holes(mask) & ~mask
|
||||
lo = np.searchsorted(values, level, side="left")
|
||||
hi = np.searchsorted(values, level, side="right")
|
||||
if hi <= lo:
|
||||
continue
|
||||
r0, r1 = int(rows_all[lo:hi].min()), int(rows_all[lo:hi].max()) + 1
|
||||
c0, c1 = int(cols_all[lo:hi].min()), int(cols_all[lo:hi].max()) + 1
|
||||
box_mask = np.pad(burned[r0:r1, c0:c1] == level, 1)
|
||||
interior = binary_fill_holes(box_mask) & ~box_mask
|
||||
if not interior.any():
|
||||
continue
|
||||
components, count = label(interior)
|
||||
for component_id in range(1, count + 1):
|
||||
whole = components == component_id
|
||||
if (constrained & whole).any():
|
||||
for component_id, obj in enumerate(find_objects(components), start=1):
|
||||
whole = components[obj] == component_id
|
||||
top, left = obj[0].start - 1 + r0, obj[1].start - 1 + c0
|
||||
bottom, right = top + whole.shape[0], left + whole.shape[1]
|
||||
if (constrained[top:bottom, left:right] & whole).any():
|
||||
continue # 안에 다른 제약이 있으면 마루가 아니다(보통의 감싸는 링)
|
||||
win = _window(whole, pad)
|
||||
component = whole[win]
|
||||
line = mask[win]
|
||||
win = (
|
||||
slice(max(top - pad, 0), bottom + pad),
|
||||
slice(max(left - pad, 0), right + pad),
|
||||
)
|
||||
line = burned[win] == level
|
||||
component = np.zeros(line.shape, dtype=bool)
|
||||
component[
|
||||
top - win[0].start : bottom - win[0].start,
|
||||
left - win[1].start : right - win[1].start,
|
||||
] = whole
|
||||
local = surface[win]
|
||||
# 그 묶음을 감싼 고리 칸 → 고리 바로 바깥
|
||||
loop = binary_dilation(component, structure=eight) & line
|
||||
@@ -164,7 +184,7 @@ def build_contour_surface(
|
||||
if len(present) < 2:
|
||||
return None
|
||||
interval_m = float(np.diff(np.array(present)).min())
|
||||
surface = build_laplace(spec, burned, None, spec.cell_m)
|
||||
surface = laplace_surface(burned, spec.cell_m)
|
||||
resolve_enclosed_interiors(burned, present, surface, spec.cell_m, interval_m)
|
||||
return surface
|
||||
|
||||
|
||||
@@ -0,0 +1,126 @@
|
||||
"""PLAN 617-2 — 공용 지형 빠른 판이 옛 길과 비트 동일인지(B04 짝).
|
||||
|
||||
① 라플라스 — B04 `build_laplace` 대 `common_util_terrain_laplace.laplace_surface` ·
|
||||
짝 · 홀 크기 · 고정 칸 많음 · 적음 · NaN 칸.
|
||||
② 폐합 채움 — 수위 상자 · 묶음 창 판 대 통격자 판(옛 셈 그대로 본뜸) · 고리가 격자 끝에 닿음 ·
|
||||
겹 고리 · 안에 다른 수위.
|
||||
|
||||
./venv/Scripts/python.exe -m pytest resources/tester/test_617_2_laplace_mirror.py -q
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import pytest
|
||||
|
||||
from B04_PreProcess.B04_PreProcess_Engine_SheetMethods import build_laplace, relax_laplace
|
||||
from B04_PreProcess.B04_PreProcess_Engine_Watershed_Grid import GridSpec
|
||||
from common_util import common_util_terrain_surface as surf
|
||||
from common_util.common_util_terrain_laplace import laplace_surface, relax_checker
|
||||
|
||||
|
||||
def _burned(shape, fixed_share, seed):
|
||||
rng = np.random.default_rng(seed)
|
||||
burned = np.full(shape, np.nan, dtype=np.float32)
|
||||
pick = rng.random(shape) < fixed_share
|
||||
burned[pick] = rng.integers(0, 12, int(pick.sum())).astype(np.float32) * 5.0 + 100.0
|
||||
return burned
|
||||
|
||||
|
||||
@pytest.mark.parametrize("shape", [(40, 52), (41, 53), (40, 53), (41, 52), (1, 7), (2, 3)])
|
||||
@pytest.mark.parametrize("fixed_share", [0.02, 0.4])
|
||||
def test_laplace_matches_b04(shape, fixed_share):
|
||||
burned = _burned(shape, fixed_share, seed=shape[0] * 100 + shape[1])
|
||||
burned.flat[0], burned.flat[-1] = 100.0, 155.0 # 수위 둘 이상
|
||||
spec = GridSpec(x_min=0.0, y_max=shape[0] * 2.0, cell_m=2.0, n_rows=shape[0], n_cols=shape[1])
|
||||
old = build_laplace(spec, burned, None, 2.0)
|
||||
new = laplace_surface(burned, 2.0)
|
||||
assert new.dtype == old.dtype == np.float32
|
||||
assert np.array_equal(old, new)
|
||||
|
||||
|
||||
@pytest.mark.parametrize("shape", [(30, 30), (31, 29)])
|
||||
def test_relax_keeps_nan_and_fixed(shape):
|
||||
rng = np.random.default_rng(3)
|
||||
base = rng.random(shape).astype(np.float32) * 50.0
|
||||
base[rng.random(shape) < 0.05] = np.nan
|
||||
fixed = rng.random(shape) < 0.1
|
||||
old, new = base.copy(), base.copy()
|
||||
relax_laplace(old, fixed, 37)
|
||||
relax_checker(new, fixed, 37)
|
||||
assert np.array_equal(old, new, equal_nan=True)
|
||||
|
||||
|
||||
def _full_grid_fill(burned, present, surface, cell_m, interval_m):
|
||||
"""옛 셈(PLAN 607-1 판) — 수위마다 통격자 채움 · 묶음마다 통격자 `components == id`."""
|
||||
from scipy.ndimage import binary_dilation, binary_fill_holes, distance_transform_edt, label
|
||||
|
||||
band_cells, constrained = surf._BAND_CELLS, np.isfinite(burned)
|
||||
limit, pad, eight = max(interval_m - 0.5, 0.5), surf._BAND_CELLS + 2, np.ones((3, 3), bool)
|
||||
resolved = 0
|
||||
for level in present:
|
||||
mask = burned == level
|
||||
interior = binary_fill_holes(mask) & ~mask
|
||||
components, count = label(interior)
|
||||
for component_id in range(1, count + 1):
|
||||
whole = components == component_id
|
||||
if (constrained & whole).any():
|
||||
continue
|
||||
rows, cols = np.flatnonzero(whole.any(axis=1)), np.flatnonzero(whole.any(axis=0))
|
||||
win = (
|
||||
slice(max(rows[0] - pad, 0), rows[-1] + pad + 1),
|
||||
slice(max(cols[0] - pad, 0), cols[-1] + pad + 1),
|
||||
)
|
||||
component, line, local = whole[win], mask[win], surface[win]
|
||||
loop = binary_dilation(component, structure=eight) & line
|
||||
ring = binary_dilation(loop, structure=eight) & ~component & ~line
|
||||
ring &= np.isfinite(local)
|
||||
if not ring.any():
|
||||
continue
|
||||
direction = 1.0 if float(local[ring].mean()) < level else -1.0
|
||||
inner = distance_transform_edt(component, sampling=cell_m)
|
||||
outer = distance_transform_edt(~(loop | component), sampling=cell_m)
|
||||
band = (outer > 0) & (outer <= band_cells * cell_m)
|
||||
band &= np.isfinite(local) & ~component & ~line
|
||||
slope = (
|
||||
float(np.mean(np.abs(level - local[band].astype(np.float64)) / outer[band]))
|
||||
if band.any()
|
||||
else 0.0
|
||||
)
|
||||
if slope > 1e-3:
|
||||
offset = np.minimum(slope * inner[component], limit)
|
||||
else:
|
||||
peak = float(inner.max())
|
||||
offset = (limit / 2.0) * (inner[component] / peak) if peak > 0 else 0.0
|
||||
local[component] = level + direction * offset
|
||||
resolved += 1
|
||||
return resolved
|
||||
|
||||
|
||||
def _ellipse(cx, cy, rx, ry, n=64):
|
||||
t = np.linspace(0.0, 2 * np.pi, n + 1)
|
||||
return np.column_stack([cx + rx * np.cos(t), cy + ry * np.sin(t)])
|
||||
|
||||
|
||||
@pytest.mark.parametrize("n_rows,n_cols", [(150, 200), (151, 199)])
|
||||
def test_window_fill_matches_full_grid(n_rows, n_cols):
|
||||
lines = [
|
||||
(_ellipse(100, 150, 80, 60), 95.0),
|
||||
(_ellipse(100, 150, 40, 30), 100.0), # 봉우리
|
||||
(_ellipse(100, 150, 12, 9), 105.0), # 겹 봉우리
|
||||
(_ellipse(300, 150, 60, 50), 105.0),
|
||||
(_ellipse(300, 150, 25, 20), 100.0), # 웅덩이(같은 높이 다른 고리)
|
||||
(_ellipse(4, 290, 30, 20), 100.0), # 격자 끝 · 모서리에 걸린 고리
|
||||
(_ellipse(396, 6, 25, 25), 110.0),
|
||||
(_ellipse(200, 60, 30, 20), 100.0),
|
||||
(_ellipse(200, 60, 8, 6), 90.0), # 안에 다른 수위 → 건너뜀
|
||||
(np.array([[0.0, 240.0], [400.0, 200.0]]), 95.0),
|
||||
]
|
||||
spec = GridSpec(x_min=0.0, y_max=n_rows * 2.0, cell_m=2.0, n_rows=n_rows, n_cols=n_cols)
|
||||
a, b, z = surf.segments_from_lines(lines)
|
||||
burned = surf.burn_segments(spec, a, b, z)
|
||||
present = sorted(np.unique(burned[np.isfinite(burned)]).tolist())
|
||||
base = laplace_surface(burned, 2.0, iterations=30)
|
||||
old, new = base.copy(), base.copy()
|
||||
count_old = _full_grid_fill(burned, present, old, 2.0, 5.0)
|
||||
count_new = surf.resolve_enclosed_interiors(burned, present, new, 2.0, 5.0)
|
||||
assert count_old == count_new >= 3
|
||||
assert np.array_equal(old, new)
|
||||
Reference in New Issue
Block a user